**Simula SpringerBriefs on Computing** Reports on Computational Physiology **13**

**Kimberly J. McCabe** *Editor*

# **Computational Physiology** Simula Summer School 2022 − Student Reports

### **Simula SpringerBriefs on Computing**

### Reports on Computational Physiology

Volume 13

### **Series Editors**

Kimberly J. McCabe Simula Research Laboratory, Oslo, Norway Andrew D. McCulloch, Institute for Engineering in Medicine, University of California San Diego, La Jolla, California, USA; Department of Bioengineering, University of California San Diego, La Jolla, California, USA

### **Managing Editor**

Jennifer Hazen, Simula Research Laboratory, Oslo, Norway

### **Editor-in-Chief**

Aslak Tveito, Simula Research Laboratory, Oslo, Norway

### **About this series**

In 2016, Springer and Simula launched an Open Access series called the Simula SpringerBriefs on Computing. This series aims to provide concise introductions to the research areas in which Simula specializes: scientific computing, software engineering, communication systems, machine learning and cybersecurity. These books are written for graduate students, researchers, professionals and others who are keenly interested in the science of computing, and each volume presents a compact, state-of-the-art disciplinary overview and raises essential critical questions in the field.

Simula's focus on computational physiology has grown considerably over the last decade, with the development of multi-scale mathematical models of excitable tissues (brain and heart) that are becoming increasingly complex and accurate. This sub-series represents a new branch of the SimulaSpringer Briefs that is specifically focused on computational physiology. Each volume in this series will explore multiple physiological questions and the models developed to address them. Each of the questions will, in turn, be packaged into a short report format that provides a succinct summary of the findings and, whenever possible, the software used will be made publicly available.

Kimberly J. McCabe Editor

## Computational Physiology

Simula Summer School 2022 − Student Reports

*Editor*  Kimberly J. McCabe Simula Research Laboratory Oslo, Norway

ISSN 2512-1677 ISSN 2512-1685 (electronic) Simula SpringerBriefs on Computing ISBN 978-3-031-25373-7 ISBN 978-3-031-25374-4 (eBook) https://doi.org/10.1007/978-3-031-25374-4 ISSN 2730-7735 ISSN 2730-7743 (electronic) Reports on Computational Physiology

Mathematics Subject Classification (2020): 92-10, 92C99

© The Editor(s) (if applicable) and The Author(s) 2023. This book is an open access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

### **Preface**

Since 2014, we have organized an annual summer school in computational physiology. The school starts in June each year and the graduate students spend two weeks in Oslo learning the principles underlying mathematical models commonly used in studying the heart and the brain. At the end of their stay in Oslo, the students are assigned a research project to work on over the summer. In August the students travel to the University of California, San Diego to present their fndings. Each year, we have been duly impressed by the students' progress and we have often seen that the results contain the rudiments of a scientifc paper.

Starting in the 2021 edition of the summer school, we have taken the course one step further and aim to conclude every project with a scientifc report that passes rigorous peer review as a publication in this new series called *Simula SpringerBriefs on Computing – reports on computational physiology*.

One advantage of this course adjustment is that we have the opportunity to introduce students to scientifc writing. To ensure the students get the best introduction in the shortest amount of time, we have commissioned a professional introduction to science writing by Nature. The students participate in a 2-day *Nature Masterclasses* workshop, led by two editors from Nature journals, in order to strengthen skills in high quality scientifc writing and publishing. The workshop is tailored to publications in the feld of computational physiology and allows students to gather real-time feedback on their reports.

We would like to emphasise that the contributions in this series are brief reports based on the intensive research projects assigned during the summer school. Each report addresses a specifc problem of importance in physiology and presents a succinct summary of the fndings (8-15 pages). We do not require that results represent new scientifc results; rather, they can reproduce or supplement earlier computational studies or experimental fndings. The physiological question under consideration should be clearly formulated, the mathematical models should be defned in a manner readable by others at the same level of expertise, and the software used should, if possible, be made publicly available. All reports in this series are subjected to peer-review by the other students and supervisors in the program.

We would like to express our gratitude for the very fruitful collaboration with Springer -Nature and in particular with Dr. Martin Peters, the Executive Editor for Mathematics, Computational Science and Engineering.

The editors of *Simula SpringerBriefs on Computing – reports on computational physiology*:

Oslo, Norway December 2022

*Kimberly J McCabe Andrew D McCulloch Aslak Tveito Jennifer L Hazen*

### **Acknowledgements**

The Simula Summer School in Computational Physiology is a team efort, with many scientists contributing their time to give lectures and advise projects for the students. We would like to thank the lecturers and project advisors for their expertise and willing participation in the course: Dr. Hermenegild Arevalo, Dr. Gabriel Balaban, Dr. Jonas van den Brink, Dr. Francisco Contijoch, Dr. Cecile Daversin-Catty, Dr. Andrew ´ Edwards, Dr. Ada Johanne Ellingsrud, Dr. Nickolas Forsch, Dr. Ingeborg Gjerde, Dr. Ilse van Herck, Dr. Glenn Lines, Dr. Molly Maleckar, Dr. Kimberly McCabe, Dr. Andrew McCulloch, Lena Myklebust, Oscar Odeigah, Dr. Marie Rognes, Dr. Joakim Sundnes, Dr. Aslak Tveito, Julie Uv, Dr. Daniela Valdez-Jasso, Dr. Samuel Wall, and Kei Yamamoto.

Administrative support for the school was provided by Dr. Kimberly McCabe, Dr. Nickolas Forsch, Dr. Rachel Thomas, Dr. Stian Engen, and Hanie Tampus. At UCSD, we received vital administrative support from Jennifer Stowe and Marlene San Pedro. Formatting support for this volume was provided by Dr. Jennifer Hazen.

The Simula Summer School in Computational Physiology is supported through the Simula-UiO-UCSD Research PhD Training Programme (SUURPh), an endeavour funded by the Norwegian Ministry of Education and Research. Additional fnancial support is derived from SIMENT, an INTPART mobility grant from the Norwegian Research Council. The school also received funding from Digital Life Norway (DLN) to support student mobility.

### **Contents**






### **Chapter 1 An Automated Cardiac Constitutive Modelling Framework with Evolutionary Strain Energy Functions**

Kristin Ludwicki\*, Leto L Riebel\*, Sophia Ohnemus\*, Frida M E Westby, Nickolas Forsch, and Gabriel Balaban \* These authors contributed equally to this work

**Abstract** Heart disease is the leading cause of mortality worldwide. Many cardiac diseases are associated with altered elastic energy relations of the heart tissue. However, the strain energy functions describing these characteristics are limited since they were designed manually for highly specifc experimental setups. In this study, we develop CHESRA (Cardiac Hyperelastic Evolutionary Symbolic Regression Algorithm), an automated constitutive modelling framework, to derive cardiac elastic strain energy functions directly from experimental data. Our results indicate that CHESRA fnds functions that reproduce mechanical tissue properties from experimental data whilst controlling function complexity. Our novel approach has the potential to fnd strain-energy functions that ft to various experimental data sets and may contribute to automatically building mathematical models to understand clinical observations of heart diseases.

Sophia Ohnemus Institute for Experimental Cardiovascular Medicine, University Heart Center Freiburg - Bad Krozingen, Medical Center and Faculty of Medicine, University of Freiburg, DE

Frida M E Westby Department of Informatics, University of Oslo, NO, Department of Computational Physiology, Simula Research Laboratory, NO

Nickolas Forsch Department of Computational Physiology, Simula Research Laboratory, NO

Gabriel Balaban (corresponding author) Department of Computational Physiology, Simula Research Laboratory, NO e-mail: gabrib@ simula.no

Kristin Ludwicki

Department of Physiology and Biophysics,Weill Cornell Medicine, New York, USA

Leto L Riebel Department of Computer Science, University of Oxford, UK

### **1.1 Introduction**

Cardiovascular disease is the leading cause of mortality worldwide, accounting for around 19 millions deaths in 2020 [1]. In many cardiac diseases, the mechanical and structural properties of the heart are altered. Examples include hypertrophic cardiomyopathy, which is characterized by a thickening of the heart walls [2] and dilated cardiomyopathy, a disease which causes the ventricular walls to thin and stretch [3]. Another example is myocardial infarction, commonly known as a heart attack, during which oxygen and nutrient supply to the cardiac tissue is interrupted [4]. Myocardial infarction results in altered mechanical loading conditions [5], structural changes, increased tissue stifness, and ultimately ventricular dysfunction [6]. A better understanding of these pathological mechanisms may improve clinical diagnosis and treatment options.

To characterize the elastic properties of the heart, stress-strain relations of the myocardium have been extensively studied experimentally. Assessment of passive shear properties of pig ventricular myocardium by Dokos et al. [7] indicates that the myocardium presents as an orthotropic material since the shear response difers along three mutually orthogonal directions. Demer and Yin [8] and Yin et al. [9] measured the stress-strain relationship of passive non-contracting canine myocardium for simultaneous biaxial stretching and observed highly nonlinear, anisotropic, and viscoelastic behavior. Novak et al. [10] performed more biaxial experiments in canine myocardium from diferent regions of the left ventricle and found that qualitatively the mechanical properties were similar, but quantitatively they difered. For example, samples from sub-endocardium and sub-epicardium of the left ventricular free wall tended to be more rigid than samples from the mid-myocardium.

In order to mathematically represent the mechanical properties of the heart tissue, a number of strain-energy functions (SEF) have been proposed. These SEF describe the potential energy density of the heart tissue depending on its deformation and have been developed based on a few key considerations. Considerations include structural metrics and other experimental results, as well as assumptions on physical properties of the system. Although the myocardium appears to be a viscoelastic material, the relaxation time of the viscoelastic response is long compared to the cardiac cycle [11]. Therefore, most SEF, with exceptions as presented in [12], model it as hyperelastic for simplicity reasons. Examples for SEF based on the assumption of transverse isotropy are given in [13, 14], whereas [11, 15, 16, 17] introduce orthotropic models. However, creating these SEF is labor intensive as they must be derived manually and, thus, only a limited amount of experimental data can be considered.

Evolutionary algorithms have the potential to automatically fnd functions that describe experimental data with minimal human guidance. This has been applied already to various research felds; for example, Doglioni et al. used evolutionary polynomial regression (EPR) to build a function for the relationship between air and water temperature [18]. Their method not only considered suitability of the discovered function to describe the given data, but also aimed to minimise the complexity of the formula. Montes et al. used a modifed EPR framework to create self-cleansing models for new sewer systems [19] and several diferent machine learning methods, including EPR, were used by Jamei et al. to model the density of hybrid nanofuids [20]. Furthermore, Javadi et al. established an EPR method to build constitutive equations describing the deformation of bridges and tunnels [21].

In this work we present the development of a Cardiac Hyperelastic Evolutionary Symbolic Regression Algorithm (CHESRA) that derives hyperelastic SEF for ventricular myocardium. In contrast to EPR algorithms, CHESRA is not restricted to polynomials, but searches the whole space of mathematical functions for SEF describing the given experimental data. Our results highlight the potential to automatically fnd SEF which reproduce experimental fndings of cardiac elastic properties of diferent species.

### **1.2 Methods**

In the following subsections we outline the individual components of CHESRA. First, in 1.2.1 we briefy introduce the physical principles of hyperelastic continuum mechanics that we use to defne the SEF. Then, in 1.2.2 we outline how we represent the SEF computationally as function trees. Subsequently, we explain how we ft the SEF to experimental data in 1.2.3. Lastly, we go through the individual steps of the evolutionary symbolic regression (ESR) algorithm in detail in 1.2.4.

### **1.2.1 Basic Principles of Hyperelastic Continuum Mechanics**

In order to defne our SEF, we build upon the framework developed by Holzapfel and Ogden for cardiac hyperelasticity [11]. Here, the structure of the myocardium is characterized by three orthonormal basis vectors: the fbre axis associated with the prevailing cell orientation, the sheet axis defned as the direction in the plane of the cell sheets perpendicular to the fbre direction, and the sheet-normal axis perpendicular to and .

Any point within the material can be described by a vector **X** = Í <sup>=</sup> ,, **e** , where **e** denote the basis vectors of the , , and axis. Upon deformation, the position of each point changes and the new position can be described by a vector **x**. The fundamental quantity describing such a deformation is the deformation gradient **F**, defned by = / with , = , , and where = det**F** = 1 for an incompressible material such as the myocardium. The right Cauchy-Green tensor

$$\mathbf{C} = \mathbf{F}^T \mathbf{F},\tag{1.1}$$

and the Green-Lagrange strain tensor

$$\mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{I}),\tag{1.2}$$

are associated with **F**, with **I** being the identity tensor.

Since these quantities do not only account for deformation but also for rotation and translation, in the following, we consider the invariants of **C** which are defned as

$$\begin{aligned} I\_1 &= \text{tr}\mathbf{C}, & I\_2 &= \frac{1}{2} \left[ I\_1^2 - \text{tr}\mathbf{C} \right] & & I\_3 &= \text{det}\mathbf{C}, \\ I\_{4i} &= \mathbf{e}\_i^T (\mathbf{C}\mathbf{e}\_i), & I\_{5i} &= \mathbf{e}\_i^T (\mathbf{C}^2\mathbf{e}\_i), & & I\_{8if} = \mathbf{e}\_i^T (\mathbf{C}\mathbf{e}\_j), \end{aligned} \tag{1.3}$$

with ≠ ∈ { , , }. Defning a SEF based on these invariants ensures objectivity. The SEF proposed by Holzapfel and Ogden [11] is given by

$$\psi = \frac{a}{2b} \exp\{b(I\_1 - 3)\} + \sum\_{i=f,s} \frac{a\_i}{2b\_i} \left\{ \exp\{b\_i(I\_{4i} - 1)^2\} - 1 \right\}$$

$$+ \frac{a\_{fs}}{2b\_{fs}} \left\{ \exp\{b\_{fs}I\_{8fs}^2\} - 1 \right\}. \tag{1.4}$$

Here, , , , , , , , and are material parameters chosen to either ft the corresponding Cauchy stress tensor

$$\boldsymbol{\sigma} = \boldsymbol{J}^{-1} \mathbf{F} \sum\_{k} \frac{\partial \boldsymbol{\psi}}{\partial I\_{k}} \frac{\partial I\_{k}}{\partial \mathbf{F}} \tag{1.5}$$

to the shear data from Dokos et al. [7], or the associated second Piola-Kirchhof stress tensor

$$\mathbf{S} = J \mathbf{F}^{-1} \sigma \mathbf{F}^{-T} \tag{1.6}$$

to the biaxial stretch data by Yin et al. [9] (see Section 1.2.3.1 for details about these experimental data sets).

### **1.2.2 Representation of the Strain-Energy Functions as Function Trees**

In CHESRA, SEF are represented as a list of nodes composing function trees. Details on function tree implementation and use cases can be drawn from [22, 23]. In our implementation, each node in the function tree has a type, a value, up to one parent, and up to two children (one left and one right). Only the root node of the function tree has no parent, further explained in 1.2.4.1. Nodes with no children are called 'leaves'. Sub-trees within the function tree represent individual terms of the SEF.

Nodes are either a symbol, an operand, or a pre-operand. Nodes of type symbol can either be a constant material parameter or an invariant. Material parameters are represented by the symbols {1, 2, 3, 4, 5, 6, 7, 8}. These material parameters are later replaced by a constant real number to ft the SEF to experimental data, as described in 1.2.3. Similarly, invariant nodes are assigned a symbol from the set of invariants {1, 2, 3, <sup>4</sup> , 4, 4, <sup>5</sup> , 5, 5, <sup>8</sup> , <sup>8</sup> , 8 }, as defned in Equation 1.3. Operand nodes can either have the value + (addition), − (subtraction), ÷ (division) or ∗ (multiplication); and pre-operands are either (to the power of 2), (euler's number '' to the power of) or - (negation). While operands require two terms, for example + , a pre-operand requires only one, for instance , which is written before the term in many coding languages, such as in python: (,2). An example function tree is shown in Figure 1.1.

Fig. 1.1: Example SEF and corresponding function tree with node types indicated. The function tree is converted to a real function by applying an in-order traversal.

If a node of type symbol, *i.e.* either a material parameter or an invariant, has a parent node, it has to be either an operand or a pre-operand to be a valid function. Similarly, symbol nodes can only have children of type operand. Pre-operands are applied to terms, and hence can only have exactly one child, which has to be of type symbol or pre-operand. Their parent has to be either an operand or another pre-operand. Operands have a parent and exactly one child, both of which cannot be another operand. If an operand is a right child, its child needs to also be a right child and vice versa, to avoid two symbols next to each other. These rules are summarised in Table 1.1. We apply brackets only around the terms pre-operators are applied to, all other operations are carried out according to their priority, e.g. ∗ before +.

Using the rules summarized in Table 1.1, to convert our function trees into functions we apply an in-order traversal, where any node is read out as: value of its left child, its own value, value of its right child. Figure 1.1 shows an example of a function and one of its possible corresponding function trees, depending on the order the terms are created in. Note that since our algorithm does not perform any searches within the tree, there is no need for it to be balanced.

Representing SEF as function trees has several benefts. Firstly, individual terms can be identifed as sub-trees, allowing for pre-operands to be applied to parts of the equation only. This would be more difcult, for example, in a linked list structure. Secondly, terms can be manipulated, swapped, and deleted easily during mutation and mating, which are explained in 1.2.4.2.


Table 1.1: Summary of node types, possible node values, and node connection rules within the function tree.

### **1.2.3 Fitting the Strain-Energy Functions to Experimental Data**

For each SEF we ft the material parameters to experimental data sets in a separate step, in order to ensure that only the material parameters difer for diferent experimental setups while the general form of the SEF is the same. In this work we consider the shear data of Dokos et al. [7] and the biaxial stretch data by Yin et al. [9], which we briefy describe in the following.

#### **1.2.3.1 Experimental Data under Consideration**

Dokos et al. [7] measured the shear stress, corresponding to the Cauchy stress tensor components versus the amount of shear in a cube of pig left ventricular myocardium sheared in the , , or plane (see Figure 1.2a). They considered diferent shear modes ( ) referring to shear in the plane in direction with ≠ ∈ { , , }.

The biaxial stretch data from Yin et al. [9] was collected in canine left ventricular tissue which was stretched simultaneously in the fbre and sheets direction. Figure 1.2b shows the measured second Piola-Kirchhof stress versus strain and respectively versus for three diferent constant strain ratios = /.

For both data sets, we used the digitizer software WebPlotDigitizer [24] to digitize the data from the graphs shown in [7] and [9].

#### **1.2.3.2 Fitting to the Experimental Data**

In order to ft the material parameters of a SEF to the shear data set we calculated the predicted Cauchy stress tensor components according to Equation 1.5. Then, we used the Python package lmft [25] to minimize the residuals,

$$res\_k = s\_{if}(\gamma\_k) - \sigma\_{if}(\gamma\_k),\tag{1.7}$$

Fig. 1.2: Experimental data sets under consideration. (a) Plot of the digitized data collected by Dokos et al. [7] in pig left ventricular myocardium for shear stress versus the amount of shear for shear ( ) in the plane in direction with ≠ ∈ { , , }. Results for ( ) and () shear are identical. (b) Visualization of the data by Yin et al. [9] for biaxial loading in the plane of canine left ventricular myocardium. The left plot shows the stress versus strain in the fbre direction and the right plot versus in the sheets direction for three diferent strain ratios = / = 2.05 (triangles), 1.02 (squares) and 0.48 (circles). In both cases we used WebPlotDigitizer [24] to digitize the data shown in [7] and [9].

between experimental data points and the corresponding calculated stress tensor component for all shear modes ( ) and all amounts of shear in the data set simultaneously.

Analogously, for ftting the material parameters to the biaxial data set we calculated the diagonal components of the second Piola-Kirchhof stress tensor according to Equation 1.6. We used the same Python package for simultaneously minimizing the residuals

$$res\_k = s\_{i\bar{i}}(E\_{i\bar{i},k}) - S\_{i\bar{i}}(E\_{i\bar{i},k}),\tag{1.8}$$

with = , , for all Green-Lagrange strain tensor components , and strain ratios = / in the data set.

### **1.2.4 Evolutionary Symbolic Regression Algorithm**

To develop CHESRA, we build upon the traditional ESR framework. First, populations of random SEF are created (see Section 1.2.4.1). Then, every equation within each iteration, or generation, is assigned an error-score based on a user-designed ftness function. Our ftness function assesses the complexity of the SEF and how well it fts to the shear and/or biaxial data set (see Section 1.2.4.3). The individuals with the best ftness are then selected, mated to generate a new population, and mutated to ensure diversity for the next generation (see Section 1.2.4.2). This multi-step process continues until the maximum number of generations is reached. A schematic representation of CHESRA is illustrated in Figure 1.3.

Fig. 1.3: The main workfow of CHESRA.

### **1.2.4.1 Setup and Initialisation**

To initialise SEF and their corresponding function trees, a maximum function length, set of invariants, and material parameters are defned. First, our algorithm creates a root node, which is always an invariant or a material parameter. Second, it picks a random length within the given maximum and extends the function tree until this chosen length is reached. The function can either be extended by adding a pre-operand in front of a term, or by attaching an operand and a new symbol (i.e., an invariant or material parameter) to an existing symbol. For the latter, the selected existing symbol may have at most one child. Details concerning function development can be found in Section 1.2.2.

#### **1.2.4.2 Mutation and Mating of Strain-Energy Functions**

We mutate SEF by randomly selecting a node from the corresponding function tree to be mutated. The selected node is then replaced by a randomly chosen element of the same type (see 1.4).

To mate two parent SEF and create two new children, our algorithm creates a deep copy of each parent function tree, chooses one random sub-tree per copy, and then swaps them over. To make sure the resulting children SEF are valid functions, we restrict the chosen sub-trees to be proper terms and hence start with a symbol or pre-operand. Alternatively, one may also choose to ensure the sub-trees to be swapped to start with a node of the same type. Our chosen sub-trees may be of diferent size and hence children SEF may be of diferent lengths than their parents (see 1.5).

Mutation and mating are important to increase population diversity in CHESRA. However, too much diversity could also prohibit convergence towards an optimal solution. Therefore, we set rates that will dictate the probability of a given individual to mate, mate, and/or mutate, mutate. These rates were set to mate = 20% and mutate = 80% in all CHESRA experiments.

Fig. 1.4: An examplary function tree representing a SEF is mutated by choosing a random item and replacing it with a random value of the same type, as highlighted by the orange boxes.

#### **1.2.4.3 Scoring Fitness of Strain-Energy Functions**

We designed a ftness function to ensure that CHESRA can penalize against poor ftting SEF with high complexity. Therefore, each SEF, or individual, in each generation is evaluated based on its complexity and how well it reproduces the results of a specifc experimental data set. It is this evaluation that allows the best individuals to be selected for mating, mutation, and population of the following generation.

Fig. 1.5: Mating between two parent SEF to generate two children SEF is implemented by swapping two random subtrees from the corresponding function trees as highlighted in the orange boxes.

To evaluate the ftness of each equation, frst we quantify its complexity as the sum of the function length eq (i.e., the number of nodes it consists of), and the number of pre-operators pre-op that it contains. The pre-operator count allows for penalization against nestedness, or the creation of functions with pre-operators within themselves.

$$Complexity = l\_{\text{eq}} + n\_{\text{pre-op}} \,. \tag{1.9}$$

Second, we assess how well each individual SEF reproduces experimental data. For this, the material parameters of each equation are frst ftted to the experimental data as described in 1.2.3.2. Then, for each equation the sum of squared errors (SSE) is calculated as the sum of the squared residuals between the experimental data points and the ft,

$$SSE = \sum\_{k} res\_k^2,\tag{1.10}$$

with as defned in Equation 1.7 for the shear data or in Equation 1.8 for the biaxial stretch data.

The fnal ftness score is then given by

$$Fitness = (\alpha \* Compplexity) + SSE,\tag{1.11}$$

where the hyperparameter ensures appropriate balance between ftting and complexity. Numerical experiments conducted to fnd a suitable can be found in Section 1.3.1.

Where CHESRA is given data sets, the ftness function can be extended to

$$Fitness = (\alpha \* Compplexity) + \sum\_{N} SSE\_N. \tag{1.12}$$

Hence, in the following section 1.3.1 we used Equation 1.11 in order to fnd a SEF that can reproduce the shear data alone, while in section 1.3.2 we used Equation 1.12 to consider both the shear and the biaxial data set.

### **1.3 Results**

### **1.3.1 Trial 1: SEF Describing the Shear Data**

In our frst trial, we ran CHESRA while ftting the SEF to the shear data set from Dokos et al. [7]. We chose the hyperparameter of the ftness function (see Equation 1.10) by running CHESRA with = 10−<sup>5</sup> , 10−<sup>3</sup> , 10−<sup>1</sup> , 1, 2, and 5 for 50 generations and 100 individuals each. 1.6 indicates that a value of = 0.1 lead to the best result as it allowed for SEF with low SSE and minimal complexity. Using this hyperparameter, the ftness score of the best SEF gradually decreased over the generations until a plateau was reached (see 1.7).

In this run, the SEF with the lowest ftness score derived using CHESRA reproduces the shear data well (see 1.8a) and is given by

$$\psi\_1 = \exp\left(p\_1 I\_{\mathbb{S}s} + p\_2 I\_{\mathbb{S}f} + 2I\_2 - I\_{4n} - p\_3\right). \tag{1.13}$$

In order to assess whether CHESRA is reproducible, we repeated this trial twice, with the best SEF given by

$$
\psi\_2 = -\left(I\_{\mathbb{S}f} - p\_1\right)^2 + p\_2 I\_{\mathbb{S}f} \left(\left(I\_{\mathbb{S}} + p\_3\right)^2 + I\_{4s}\right)^2,\tag{1.14}
$$

$$\psi\_3 = \exp\left(p\_1^2 I\_{\xi f} + (p\_2 + p\_3 - 1)I\_{8fs} + p\_4 I\_{4f} + I\_1 + 2I\_2 - I\_{4n} + p\_S\right). \tag{1.15}$$

The equations derived from the three trial runs have unique material parameters that can be found in 1.2. 1.8 shows a good ft to the shear data for all three SEF found.


Table 1.2: Material parameters to ft the SEF defned in Equation 1.13 (trial run 1), 1.14 (trial run 2), 1.15 (trial run 3) to the shear data set by Dokos et al. [7]

Fig. 1.6: SSE vs. complexity for diferent hyperparameters in the ftness function (see Equation 1.11). Six runs of CHESRA were completed in which each was run with one of the unique values listed in the legend. For experiment 1, we chose = 0.1 since it lead to SEF with low SSE and function complexity.

Fig. 1.7: Evolution of 100 individual SEF for 50 generations. Each point represents the ftness of the best SEF from each generation.

Fig. 1.8: In our numerical experiment 1, CHESRA found a SEF that reproduces the shear data of Dokos et al. [7] in all 3 trial runs. Markers represent the experimental data points whereas solid lines show the ft of the SEF derived using CHESRA.

### **1.3.2 Trial 2: SEF Describing the Shear and Biaxial Data**

As a next step, we extended CHESRA to fnd a SEF that reproduces the results of multiple experimental data sets. Therefore, we ran CHESRA with the same setup as in trial 1, but updated the ftness function to incorporate the SSE of both the shear and the biaxial data set, as described in Equation 1.12. Moreover, in this numerical experiment we used = 0 in order to assess quality of ft without penalizing for equation complexity. In addition, = 0 allowed convergence in a shorter number of generations as ftting to two diferent data sets is a challenging optimization problem.

The SEF found using CHESRA is given by

$$
\psi = \left( p\_1 p\_2 I\_{4f} - I\_{4f}^2 + \frac{p\_1}{I\_{5n}} - I\_{4s} \left[ p\_1 p\_2 I\_{5n} I\_{4f} - p\_3 p\_4^2 I\_1 + p\_5 + p\_6 + p\_7 \right]^2 \right)
$$

$$
$$

$$
\left( \begin{array}{c} \\ \end{array} \right)^2, \tag{1.16}
$$

where the material parameter values are provided in Table 1.3. Figure 1.9 shows that the equation derived in experiment 2 has a better ft to the biaxial data than to the shear data.


Table 1.3: Material parameters for ftting the SEF defned in Equation 1.16 to the shear and biaxial data sets.

Fig. 1.9: In experiment 2, CHESRA was extended to fnd a single SEF aimed to reproduce both the shear data from Dokos et al. (a) and biaxial stretch data from Yin et al. (b). Solid lines show the ft of the SEF found using CHESRA to the experimental data represented by the markers.

### **1.4 Discussion**

CHESRA automatically generates functions which reproduce experimental data of cardiac hyperelastic properties. The results of trial 1 indicate that CHESRA can derive SEF that ft the shear data from Dokos et al. [7] well. However, the quality of ft was not equal for all shear directions. Specifcally, in all trial runs the derived SEF did not ft to the shear data for the (nf) and (ns) direction, as shown in Figure 1.8. As the magnitude of shear stress is much lower in those cases, it is likely that the ftness function neglects these shear modes. Previous studies have used coefcient of determination in the ftness function [18, 21], which will be a key consideration in our future work.

To extend CHESRA to multiple experimental data sets, we included both shear [7] and biaxial data [9] in trial 2. Here, the biaxial data was ftted at the expense of the shear data. A possible cause for this could be a non-optimal individual and generation number. Since more data sets were given to the algorithm, optimization is challenging in only 50 generations. Therefore, increasing the number of generations and individuals may allow for better convergence and a more universal SEF to be found.

Ensuring that CHESRA derives simple equations was important to this study. Therefore, we included the hyperparameter to the ftness function to ensure penalization against long equations with nesting. While this is a valuable proof-of-concept, it has limitations. Specifcally, it does not take computational time into consideration. Long but simple equations might be penalized more compared to shorter, but more complex or nested functions. The time it takes for an equation to be calculated may be a better measure of simplicity and will be explored in next steps.

Moreover, the current version of CHESRA can only generate SEF using quadratic and exponential functions; however, this could be widened to a larger function space. In addition, CHESRA may be extended to even more experimental data sets. Future work may also investigate to what extent the functions generated by CHESRA fulfll the mathematical requirements of SEF, such as convexity and strong ellipticity [11]. Lastly, it would be interesting to compare our method to a Multivariate Adaptive Regression Spline (MARS) model since studies by Jamei et al. [20] suggest that this can be more accurate than ESR.

### **1.5 Conclusions**

In this work we present CHESRA, an ESR framework for automatically deriving strain-energy functions of the myocardium. Our algorithm allows for controlling function complexity and inputting multiple experimental data sources.

### **References**


regression–multi-objective genetic algorithm strategy. *Urban Water Journal*, 17(2):154–162, 2020.


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 2 Electromechanical** *In Silico* **Testing Alters Predicted Drug-Induced Risk to Develop Torsade de Pointes**

Anna Busatto, Jonathan Krauß, Evianne Kruithof, Hermenegild Arevalo, and Ilse van Herck

**Abstract** Torsade de Pointes (TdP) is a type of ventricular tachycardia that can occur as a side efect of several medications. The Comprehensive *in vitro* Proarrhythmia Assay (CiPA) is a novel testing paradigm that utilizes single cell electrophysiological simulations to predict TdP risk for drugs that could potentially be used clinically. However, the efects on mechanical performance and mechano-electrical feedback are neglected. Here, we demonstrate that including electromechanical simulations in CiPA testing can provide additional insights into the predicted drug-induced TdP risk. In this work, we analyzed six drugs, namely fecainide, ibutilide, metronidazole, mexiletine, quinidine and ranolazine. We compared previously classifed risks (low, intermediate, high) with our fully coupled electromechanical simulation results based upon the action potential, the electromechanical window, and the maximum active tension [1]. For ranolazine and metronidazole the predicted risk changed from low to intermediate and intermediate to high, respectively. For the latter, while electrophysiological markers indicated a low risk, the active tension decreased by 58% which can result in a loss of heart function. Therefore, adding mechanics to CiPA testing results in an altered prediction of drug-related TdP risk.

Anna Busatto

Jonathan Krauß Institute of Biomedical Engineering, Karlsruhe Institute of Technology, DE

Evianne Kruithof Eindhoven University of Technology, NL

Hermenegild Arevalo Department of Computational Physiology, Simula Research Laboratory, NO

Ilse van Herck (corresponding author)

Department of Biomedical Engineering, University of Utah, USA, Scientifc Computing and Imaging Institute, University of Utah, USA

Department of Computational Physiology, Simula Research Laboratory, NO e-mail: ilse@ simula.no

### **2.1 Introduction**

Torsade de Pointes (TdP) is a form of abnormal heart rhythm often preceded or caused by a prolonged QT-interval of the ECG [2]. Due to the high risk of unwanted and dangerous side-efects, promising new drugs with a potential TdP risk have been excluded from the research and drug-development pipeline. This, in combination with the long and complicated process of drug approval, led to a lack of new drugs on the market for cardiac disease. Therefore, the Comprehensive *in vitro* Proarrhythmia Assay (CiPA) is performed to better predict the TdP risk for drugs that could potentially be used in the clinic. However, some of the tests can result in inaccurate predictions of the actual outcomes and lead to restrictive results, limiting the number of drugs allowed on the market.

There has been ongoing research on improving these models to increase the predictability of drug-induced efects in the heart [3]. Llopis *et al.* proposed CiPA testing using a population of models approach to obtain more accurate risk predictions by simulating the efects of altered ion channel functions on the electrophysiological behavior of cardiac cells [1].

The electromechanical window () is a biomarker for TdP risk introduced by Passini *et al.* for *in silico* drug testing [4]. Clinically, is defned as the diference between the duration of electrical and mechanical systole [4]. However, for the *in silico* test, they defned the as the time diference between calcium transient duration at 90% repolarization ( 90) and action potential duration at 90% repolarization (90) [4]. This assumption makes the simulation computationally less expensive and complex, but this simplifcation reduces accuracy in the mechanical component.

In previous work, the efects of electromechanical feedback are neglected and no insight into the mechanical function is obtained. This feedback alters the electrophysiological results, especially the calcium transient, and afects the biomarkers used to predict TdP risk. However, the efect of drugs on the mechanical performance can be investigated with electromechanical simulations. Therefore, in this work, a fully coupled electromechanical simulator SimCardEMS is used to perform CiPA tests [5]. Here, the cardiac mechanical function can be analyzed in terms of active tension, tissue deformations and where the latter can be altered due to electromechanical feedback. We performed several tests to compare our predicted TdP risks to previous classifcations such as the ones found in Llopis *et al.* [1].

### **2.2 Methods**

To assess TdP risk, we focused on selected biomarkers, namely maximum active tension a,max [kPa], <sup>90</sup> [ms], <sup>90</sup> [ms], and the [ms]. Using these metrics, we analyzed six diferent drugs from CiPA: fecainide, ibutilide, metronidazole, mexiletine, quinidine, and ranolazine with known cardiac efects. Each of the drugs has been classifed as high, intermediate, or low risk for TdP in Llopis *et* *al.* [1]. Here, the indicated risks are compared with the risks predicted by a fully coupled electromechanical model to investigate whether the mechanical biomarkers should also be taken into account for risk classifcation.

### **2.2.1 Drug Implementation**

The cardiac efects of the drugs are tested for three diferent plasma concentrations (). For each of the concentrations, the ion channel conductances () are multiplied with scaling factors () representing the efect of the tested drug. The values are calculated using the drug , the inhibitory concentration (50) and Hill coefcient (ℎ) value of the targeted ion channel:

$$SF = \frac{\text{g}\_{\text{drug}}}{\text{g}} = \left[1 + (\frac{PC}{IC\_{S0}})^h\right]^{-1} \tag{2.1}$$

Table 2.1: Ion channel scaling factors (SF) for quinidine using three diferent plasma concentrations (). Scaling factors for seven ion channels were calculated with Equation 2.1.


Equation 2.1 is derived from the standard Hill equation [1]. The , <sup>50</sup> and ℎ values are taken from Llopis *et al.*. The values are calculated for the considered ion channels Kr, Na, NaL, CaL, Ks, K1 and to at three diferent s which are half (0.5), normal (1.0) and double (2.0) the efective free therapeutic plasma concentration [1]. As an example, specifc values used to simulate the efect of quinidine are shown in Table 2.1.

### **2.2.2 SimCardEMS Simulations**

The calculated values are fnally used as input for the SimCardEMS solver, which is implemented in FEniCS [6], to simulate the efects of drugs on an endocardial tissue block. We solved the electrophysiological O'Hara-Rudy model which is strongly coupled with the mechanical Land model via intracellular calcium concentration and stretch activated channels [7, 8]. The active tension was scaled by a factor of fve. The Holzapfel model was used to describe the tissue material properties, where we increased the stifness by a factor of ten [9].

The simulations were performed on a 6x3x3 mm mesh representing an endocardial tissue slab with 1 mm and 0.5 mm spatial resolution for the mechanical and electrophysiological parts, respectively (Fig. 2.1). One corner of the mesh was fxed, and each of the connected planes were fxed in their respective normal directions, allowing the simulation of free contraction. The fbers were aligned in the longitudinal direction of the mesh and the tissue was activated by stimulating the entire mesh.

To reach steady state, 40 beats of 1000 ms each were simulated and used as control. Thereafter, the drugs were added into the simulation by loading the calculated s. Each of these simulations was run either for fve or ten beats (1000 ms each) depending on when steady state was reached; the goal was to investigate the efect of each tested drug on the electrophysiological and mechanical behavior of the cells. The biomarkers used for the analysis were extracted from the center of the mesh for the last beat of each simulation.

### **2.2.3 Evaluation**

The TdP risk classifcations were based upon maximum active tension a,max, 90, 90, and where

$$EM\_{\mathcal{W}} = CaTD\_{\mathcal{W}0} - APD\_{\mathcal{W}0} \tag{2.2}$$

Fig. 2.1: (Left) Mechanical mesh 6x3x3 mm with 1 mm resolution. (Right) Electrophysiological mesh 6x3x3 mm with 0.5 mm resolution. Three planes are fxed in their normal direction to allow free contraction of the material.

The metrics used for evaluation are visualized in Fig. 2.2. TdP risk was classifed as high or intermediate respectively in the following situations:


If none of these criteria were met, the drug was considered safe and classifed as low risk.

Fig. 2.2: Transmembrane voltage, calcium concentration, and active tension for control. The metrics used to evaluate the results are indicated in the fgure with dashed lines.

### **2.3 Results**

The biomarkers extracted for all simulations are shown in Table 2.2. For the control simulation we found a , of 55.0 kPa, <sup>90</sup> of 258 ms and of 134 ms. The associated risk classifcation thresholds for , were 44.0 kPa (intermediate) and 27.5 kPa (high). For 90, the thresholds were 273 ms (intermediate) and 284 ms (high), while for they were 121 ms (intermediate) and 107 ms (high).

### **2.3.1 Action Potential Duration**

The resulting action potentials for control as well as the six examined drugs at 1.<sup>0</sup> are shown in Fig. 2.3. Ibutilide, quinidine and fecainide have <sup>90</sup> values of 597 ms, 446 ms and 352 ms respectively which are above the threshold for high risk drugs. Ranolazine has an <sup>90</sup> of 283 ms at 1.<sup>0</sup> which fts in the range of intermediate risk drugs. <sup>90</sup> values for metronidazole and mexiletine are below

Table 2.2: Overview of the extracted biomarkers, a,max [kPa], <sup>90</sup> [ms], <sup>90</sup> [ms], [ms] and TdP risk classifcations for all performed simulations. 0.5, 1.<sup>0</sup> and 2.<sup>0</sup> results are given in this order for each of the six tested drugs.


the intermediate threshold and should therefore be classifed as low risk if only the <sup>90</sup> is considered.

All the chosen drugs were then examined for 0.<sup>5</sup> and 2.<sup>0</sup> as well. An example of the diferences of action potential shape for quinidine based on the PC is shown in Fig. 2.4. A concentration dependent efect is observed for <sup>90</sup> as it is the smallest for 0.<sup>5</sup> and increases with an increase in PC. An overview of the <sup>90</sup> values for each of the examined drugs at 0.5, 1.0, and 2.<sup>0</sup> is given in Table 2.2.

### **2.3.2 Electromechanical Window**

Looking at the values, high risk drugs were classifed by 20% decrease compared to control; fecainide, ibutilide, and quinidine returned values below this threshold, indicating high risk, while metronidazole and mexiletine returned values above it. Finally, for ranolazine the decrease was slightly less than 20% which, in combination with the remaining biomarkers, led us to classify it as an intermediate risk drug. For all the drugs, the TdP risk classifcation based on the was the same as <sup>90</sup> risk classifcations. In four out of the six analyzed drugs, the pro-

Fig. 2.3: Action potentials for control and six tested drugs at 1.0. Based on these results, ibutilide, quinidine and fecainide are classifed as high risk. Ranolazine is classifed as intermediate risk for this drug concentration.

gressively decreased with higher drug concentration compared to control. However, addition of metronidazole at 0.<sup>5</sup> increased the initially while for higher concentrations of metronidazole (2.0) the decreased again. The addition of mexiletine changed the by only 3 ms while both the <sup>90</sup> and <sup>90</sup> remained similar to control. All results are shown in Table 2.2.

### **2.3.3 Maximum Active Tension**

Applying metronidazole, the maximum active tension decreased from 55.0 kPa to 28.8 kPa, 22.9 kPa, and 17.8 kPa for 0.5, 1.<sup>0</sup> and 2.0, respectively. This represents a decrease in , of 48%, 58%, and 68%. For ranolazine, the , increased by 45% at 1.<sup>0</sup> compared to control. In the remaining drugs, the maximum active tension remained similar or slightly increased compared to control.

Fig. 2.4: Action potential for control and quinidine at 0.5, 1.0, and 2.0. An increase in PC highlights the drug concentration efect on the action potential duration.

### **2.3.4** *In Silico* **and** *in Vitro* **TdP Risk Classifcations**

We classifed all six drugs using the evaluation process described in section 2.2.3. Flecainide, ibutilide, and quinidine were determined to be high risk both in the previous classifcation and in our electromechanical model (based on <sup>90</sup> and values). Based on the maximum active tension, the risk of metronidazole was changed from intermediate (CiPA indication) to high risk. On the other hand, following the <sup>90</sup> and values it would have been classifed as low risk. Mexiletine was classifed as low in the CiPA as well as in our study. Lastly, ranolazine was previously predicted as low risk but was classifed as intermediate risk using our model based on both the <sup>90</sup> and values. Even though the increase in , for ranolazine was the largest of the analyzed drugs, potentially having an efect on contraction and tissue stress, increases in a,max were not considered in TdP risk classifcation.

### **2.4 Conclusions**

In this study we analyzed six drugs, fecainide, ibutilide, metronidazole, mexiletine, quinidine, and ranolazine. Three of these drugs were classifed as high risk, one as intermediate risk, and two as low risk according to CiPA. Using our fully coupled electromechanics model, four drugs were classifed as high risk, one as intermediate risk and one as low risk. An overview of this distribution is given in Table 2.3. After performing our analysis, we concluded that some drugs may belong to diferent risk categories depending on the individual parameters considered. For example, ranolazine was previously classifed as low risk; in our simulation, both the <sup>90</sup> and biomarker values for 1.<sup>0</sup> fell in the range of intermediate risk, while for 2.<sup>0</sup> it was predicted high risk. In general, for all six drugs, both the <sup>90</sup> and the biomarkers indicated the same risk categories. However, the risk classifcation based on , difered from the other biomarkers. For metronidazole, a low TdP risk was observed based on the electrophysiological biomarkers. Despite that, there was a decrease of 58% in , which can lead to severe problems in heart function.

Table 2.3: Predicted versus *in vitro* risk classifcations for the analyzed drugs. Out of the six tested drugs, two drugs were classifed diferently based on the electromechanical results.


We conclude that the addition of a mechanical biomarker, the maximum active tension, is valuable in the prediction of cardiac risk stratifcation. Additionally, using the electrophysiological biomarkers, SimCardEMS succeeded in predicting the same TdP risk category as previously determined for four of the six drugs.

2 Electromechanical Testing in TdP 29

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 3** *In silico* **Investigation of Sex-Specifc Osteoarthritis in Human Articular Chondrocytes**

Khoa Ngo, Nathaniel T Herrera, Milda Folkmanaite, Kei Yamamoto, and Mary M Maleckar

**Abstract** Osteoarthritis (OA), a progressive degenerative disease of cartilage in joints, is the most common cause of chronic disability in older adults. While OA is mostly considered an age-related pathology, women have a 1.5-fold higher risk of developing OA relative to men and experience more severe symptoms. Yet, they remain underrepresented in musculoskeletal research and clinical trials. Responsible for cartilage formation, articular chondrocytes experience physiological changes in OA, but the functional implications of such alterations remain largely unexplored due to difculties in acquiring the data experimentally. Through reparameterization, we expand a mathematical chondrocyte model to investigate sex-specifc OA pathogenesis. We performed sensitivity analysis to address the impact of ion channel activity in healthy and OA chondrocyte populations. Simulations show that in healthy female chondrocytes, the resting membrane potential is more depolarized than in healthy male chondrocytes, suggesting potential sex-specifc emergent physiological diferences in articular chondrocytes. In both sexes, the resting membrane potential of healthy chondrocytes is most sensitive to − , −, and −, but in OA it depolarizes and becomes sensitive to , and −. Developed and evaluated against experimental data, our articular chondrocyte OA electrophysiolog-

Khoa Ngo Department of Physiology and Membrane Biology, University of California, Davis, USA

Nathaniel T Herrera Department of Pharmacology, University of California, Davis, USA

Milda Folkmanaite Department of Physiology, Anatomy and Genetics, University of Oxford, UK

Kei Yamamoto Department of Computational Physiology, Simula Research Laboratory, NO

Mary M Maleckar (corresponding author)

Department of Computational Physiology, Simula Research Laboratory, NO e-mail: mmaleck@ simula.no

ical model can be used to further study OA pathology and sex-specifc pathological OA changes.

### **3.1 Introduction**

Mammalian chondrocytes are found within intervertebral discs and articular cartilage, where they maintain the extracellular matrix (ECM) and produce the cartilage matrix [1]. The ECM mainly consists of collagen fbers that transmit force and dissipate energy, proteoglycans that withstand compressional forces, and synovial fuid to reduce friction between bones. In normal function, healthy chondrocytes respond to outside stimuli and damaged articular cartilage, where they increase production of specifc ECM complexes [2]. Functional, mature, healthy articular chondrocytes are necessary for the maintenance of cartilage and therefore, physiological joint motion [3]. Chondrocytes function to maintain a homeostatic environment in the articular cartilage. Pathophysiologic conditions can lead to the development of osteoarthritis (OA); however, the mechanism by which this occurs is poorly understood. Osteoarthritis is the most common form of arthritis and is a disease characterized by the degeneration of cartilage and underlying bone [4]. Women are 1.5x more likely to develop OA, and experience more severe symptoms as compared to men [5, 6, 7]. In OA, chondrocytes secrete increased levels of infammatory cytokines, actively produce proteoglycans and collagen type II to recover the degeneration of the ECM, and become hypertrophic [8]. The present study on human, sex-specifc OA is modeled after the specifc cell physiology of the human articular chondrocyte.

Recent research suggests that maintaining a stable resting membrane potential is essential for chondrocytes in articular cartilage to be able to withstand pressure and force changes [9]. Despite chondrocytes being non-excitable cells, they undergo robust ion-channel mediated changes in response to OA. This is further supported by experiments where ion-channels contributing to the resting membrane potential were blocked, resulting in a decrease in the production of matrix mRNAs, proteins and glycosaminoglycans [10, 11]. Additionally, channel blockers inhibit chondrocyte proliferation and increase apoptosis [12, 13]. During OA, chondrocytes have reduced capability to respond to changes in the extracellular environment, specifcally with respect to ion transport; this is recognized by i.e. a defciency in volume regulation [14]. Chondrocyte diameters vary between 7 and 30 µm, making technical electrophysiological experimental studies difcult [15]. Therefore, there is limited experimental data investigating the electrophysiological diferences of chondrocytes in normal physiology and pathophysiology.

Therefore, to gain a better understanding of the human chondrocyte channelome in the context of OA, we have developed an OA chondrocyte model by further expanding a previously established mathematical model of the human articular chondrocyte [3, 16]. Mathematical models can play a key role in the dearth of physical data and in hypothesis generation, among other utilities. In the absence of constraining data, models can be used to generate valuable hypotheses about the underlying processes at work, which can then be tested in carefully designed experiments and/or against data from other sources. This initial model of the human articular chondrocyte features various ion channels describing the resting membrane potential and ion movement across the cell membrane. Previous iterations of the model focused on potassium fuxes across the cell membrane, and included functions specifcally for the integration of the sodium/potassium pump, an active transport pump that establishes an electrochemical gradient across the cell membrane. In this study, the model is expanded to include sex-specifc OA by implementing the diferences in ion channel conductances observed between males and females, as well as alterations of ion channel expression in OA disease. Populations of models using a log-normal distribution were also implemented to accurately represent biological variability, and subsequent parameter sensitivity analysis was performed to identify the currents that infuence chondrocyte resting membrane potential to the greatest degree in both health and disease.

In summary, we have developed a mathematical OA chondrocyte model to gain a better understanding of OA pathogenesis through exploratory simulations, driving forward hypothesis development in the context of relative scarcity of human chondrocyte electrophysiological data across the channelome. As OA is more prevalent in females than males, we also implemented sex diferences in our updated model. Our main goals in the present study were to: (i) update and scale the baseline male model (1) to generate a (2) disease OA-male model, (3) a control female model, and (4) a disease OA-female model, (ii) utilize a population of models approach to determine which ion channels were the main contributors to resting membrane potential in human articular chondrocytes during normal function and in OA.

### **3.2 Methods**

To study the efects of OA on the resting membrane potential, we employed the chondrocyte model and simulation protocols developed by Maleckar *et al.* [3] and further improved by Fischer-Holzhausen *et al.* [17]. The model components are shown in Figure 3.1. has been scaled to be within 0.1-0.6 of the original value, corresponding roughly to 0.9–5.35 pA/pF, as prior work has shown this range is more likely to present a physiologically-relevant model regime [17].

Following the Hodgkin-Huxley formalism, the cell membrane is modeled as a capacitor coupled in parallel with ion channels represented by resistors. The timeevolution of the transmembrane potential is given by the summation of all electrogenic transport processes.

$$C\_m \frac{dV}{dt} = -(I\_{KDR} + I\_{NaK} + I\_{NaCa} + I\_{Ca-ATP} + I\_{K-ATP})\tag{3.1}$$

$$+ I\_{K2pore} + I\_{Nab} + I\_{Kb} + I\_{Clb} + I\_{KCa})\,,$$

Fig. 3.1: Schematic illustration of the model adapted from Maleckar *et al.* [17], showing the principal ion-selective channels, exchangers, and pumps expressed in human chondrocytes.

where is the chondrocyte capacitance (6.3 pF). A set of three transport equations gives the time-derivatives of the intracellular ionic concentrations of <sup>+</sup> , + , and 2+

$$\frac{d\left[Na^{+}\right]\_{i}}{dt} = -\frac{I\_{Nab} + 3I\_{NaK} + 3I\_{NaCa} - I\_{NaH} + I\_{Clb}}{vol\_{i} \times F},\tag{3.2}$$

$$\frac{d\left[K^{+}\right]\_{i}}{dt} = -\frac{I\_{Kb} - 2I\_{NaK} + I\_{KDR} + I\_{K2pore} + I\_{KCa} + I\_{K-ATP}}{vol\_{i} \times F},\tag{3.3}$$

$$\frac{d\left[Ca^{2+}\right]\_i}{dt} = -\frac{I\_{Ca-ATP} - 2I\_{NaCl}}{2vol\_i \times F} - 0.045 \times \frac{dO\_c}{dt} \,, \tag{3.4}$$

where is the internal volume of the chondrocyte (0.005884 mL), is the Faraday constant (96,485 C/mol), and is the fraction of intracellular calmodulin bound to 2<sup>+</sup> . The coefcients in front of the exchangers come from their respective ionic exchange ratios. For example, the NaK exchangers exchange three <sup>+</sup> ions for two + ions.

### **3.2.1 Modeling Ionic Changes Induced by Osteoarthritis in Male and Female Chondrocytes**

To gain mechanistic insight into the pathogenesis of OA, we implemented changes in ion channel expressions as measured by Karlsson *et al.* [18]. We further expanded our investigation by integrating sex-specifc diferences in ion channel activity into our model. Due to the lack of experimental data on female chondrocytes, we compiled a list of sex-specifc diferences in ion channel expression based on mRNA and ion channel subunit expression data collected from epicardial cardiomyocytes. In our model, we modeled the efect of changes in ion channel expression on their conductances (), as displayed in Table 7.2.


Table 3.1: Sex-specifc subcellular ion channel conductances relative to the male baseline model

Implemented as change to the Q10 parameter which defnes the efect of temperature on ion channel kinetics [19, 20].

### **3.2.2 Generation of a Population of Models**

To capture biological variability and investigate parameter sensitivity, we used the published chondrocyte model [3, 17] as a baseline to build our population of 1,000 models by randomly modifying parameters corresponding to conductances of ion current and maximal rates of ion transports. A relatively large sample size of 1,000, typically used in similar population-based studies [21, 22], was chosen to allow even small efects, often indistinguishable in small sample sizes, to reach statistical signifcance. Prior studies have shown that biological parameters, when randomly selected, are distributed according to log-normal frequency curves instead of normal ones [23]. For each channel in the model, we varied its conductance by sampling from a log-normal distribution centered around their respective baseline value with the standard deviation set to 0.15.

### **3.2.3 Model Parameter Sensitivity Analysis**

We performed parameter sensitivity analysis on the generated model population to assess the sensitivity of the model output (e.g., the resting membrane potential) in response to channel conductance variations. The degree to which perturbing a set parameter can infuence the results are quantifed as a set of regression coefcients, which are obtained by performing multivariable partial least squares regression on the dataset. Partial least squares regression was chosen over standard multivariable regression as many independent parameters (e.g., ion channel conductances) were used to predict a smaller set of dependent variables (e.g., the resting membrane voltage) [24]. This methodology was frst used in cardiac electrophysiology by Sobie *et al.* [25] and has been widely used in other felds of biology [21, 26].

### **3.3 Results**

### **3.3.1 Modeling the Impact of Sex-Specifc OA on Resting Membrane Potential**

In Figure 3.2, we implemented the experimentally-observed changes in ion channel expression in control and OA chondrocytes, accounting for sex diferences as described in Table 7.2. Figure 3.2 shows that the healthy female chondrocyte maintains a more depolarized resting membrane potential (-58.21 ± 5.75 mV) as compared to the healthy male chondrocyte (-69.11 ± 4.71 mV). Electrophysiological diferences instigated by OA result in depolarization in both male (from -69.11 ± 4.71 mV to -53.87 ± 3.21 mV) and female (from -58.21 ± 5.75 mV to -49.03 ± 1.34 mV) chondrocytes. This overall depolarizing efect of OA remodeling on the resting membrane potential has previously been reported in experiments [27].

Fig. 3.2: Resting membrane potential in control and OA chondrocytes for (A) male and (B) female models. Changes in ion channel expression in OA induce resting membrane depolarization.

### **3.3.2 Population of Models**

Shown in Figure 3.3, and as outlined previously, we further expanded our analysis to investigate population variability by generating 1,000 models for each condition with channel conductances varied according to a log-normal distribution. For all 1,000 parameter sets, the steady-state solutions were reached within the simulation period of 50,000 seconds. Figure 3.3 shows that in male and female OA chondrocytes, intracellular sodium quickly becomes depleted and remains very small, close to zero, for the rest of the simulation. We calculated the mean resting membrane potentials for all conditions and compared them against the reported experimental values as measured in experiments, shown in Figure 3.4. In all cases, the membrane potentials of OA chondrocytes are relatively more depolarized compared to their healthy counterparts. Comparing male and female chondrocytes, although the control resting membrane potential values slightly difer (females are more depolarized compared to males), OA male and female chondrocytes depolarize to a similar value (approximately -51 mV).

### **3.3.3 Parameter Sensitivity Analysis**

Sensitivity analyses were performed on the chondrocyte populations to reveal how ionic modulations induced by OA might afect the sensitivities of underlying electrophysiological properties in the models. Shown in Figure 3.5, in the healthy male chondrocyte population, , − , −, and − are predicted to have the most critical impact on the resting membrane potential. In this case, − is the strongest depolarizing current and − is the strongest hyperpolarizing current. In male chondrocytes expressing OA, only and − remain the most infuential currents. Compared to the healthy case, the impact of − and − are severely reduced in OA, while gains a 1.6-fold depolarizing infuence.

The healthy female chondrocyte population shares some similarities with their male counterparts, as the resting membrane potential is still the most sensitive to , − , −, and −. However, − now plays a slightly more dominant role as a hyperpolarizing current. Compared to male chondrocytes, gains a more signifcant role in establishing the resting membrane potential than in males (2.5-fold increase in importance). is the strongest depolarizing current and − is the strongest hyperpolarizing current. Surprisingly, in female chondrocytes expressing OA, the majority of the currents no longer have as much of an impact on the resting membrane potential. Nonetheless, and − each still retain a minor infuence.

Fig. 3.3: Membrane voltage and internal ion concentrations for (A) male and (B) female population.

### **3.3.4 Inhibiting Restores Normal Resting Membrane Potential in OA Chondrocytes**

Per sensitivity analysis results, has the largest impact on depolarizing the resting membrane potential (i.e., having positive regression coefcients) in both male and female chondrocytes in OA. To restore the membrane potential back to near control values, we applied current block treatment then simulated the membrane potential in male and female OA populations of 1,000 models each. Figure 3.6 shows that applying 45% block to male and 55% block to female OA chondrocytes restored the mean resting membrane potential of the population back to the control level. Notably, previous clinical studies have revealed that

Fig. 3.4: Mean resting membrane potential in control and OA chondrocytes for male and female populations (1,000 models each) compared against experimental data [27]. To account for diferences in data collection methodology, experimental resting membrane potential values were scaled to match the mean of the male population.

Fig. 3.5: Parameter sensitivity analysis on (A) male and (B) female population of sex-specifc chondrocyte models. The regression coefcients measure the relative impact of a particular ionic current on the resting membrane potential. A positive regression coefcient indicates the associated current is a depolarizing current, while a negative regression coefcient indicates a hyperpolarizing current.

blockers Ouabain and Digoxin protected against OA and relieved OA-associated pain [28].

Fig. 3.6: Resting membrane potential in control, OA, and OA chondrocytes with partial block of for male and female model population (1,000 models each). Blocking (45% for male and 55% for female population) helped restore the resting membrane potential back to normal levels in OA chondrocytes.

### **3.4 Discussion**

OA, a debilitating disorder afecting 1 in 7 individuals in the United States in their lifetime [29], has no current curative treatment. The prevalence, as well as disease severity, is more pronounced in female patients [30]. It is known that chondrocytes play an essential role in maintaining healthy cartilage, and that their function is impaired in disease [31]. As experimental investigation of human chondrocytes is challenging due to the lack of human control samples, we have further expanded a previously-developed mathematical model of a chondrocyte channelome by introducing experimental changes observed in ion channel electrophysiology, as well as temperature and capacitance [3]. The updated model is user-friendly and freely accessible (github.com/k-ngo/Chondrocyte).

Simulations revealed that the chondrocyte membrane becomes depolarized in OA, which agrees with experimental data [32]. This suggests that the changes in ion channel expression together with temperature and capacitance are sufcient to capture the main pathological alterations observed in OA chondrocytes related to their membrane potential. In the future, it remains important to further validate the model experimentally, as well as to investigate the primary cause of the alterations in ion channel expression observed in OA chondrocytes. These are likely to involve infammatory processes associated with OA [33].

Comparisons between male and female chondrocyte models reveal intriguing phenotypic diferences. First, the resting membrane potential is more depolarized in females compared to males in controls. The main diferences between male and female healthy chondrocyte models include reduced (by 25%), − (by 95.2%), (by 2.8%) and (by 5.7%). The investigation of the main factors important for maintaining membrane potential visualized in Figure 3.5 revealed that, in the present model, calcium ATPase and background sodium and potassium currents contribute the most for maintaining healthy resting membrane potential in males, with an addition of in females. It could be that the relative depolarization of the female chondrocyte membrane driven by the aforementioned changes in the model predisposes female patients to OA, although additional experimental data is needed to further investigate this hypothesis.

OA chondrocyte model simulations show that resting membrane potential is depolarized in disease, following experimental data [32]. This efect is seen in both male and female OA models. Our results also suggest that chondrocyte resting membrane potential in female patients could be more depolarized compared to that in males, which could lead to higher depolarization in OA (mean resting membrane potential is -49.03 ± 1.34 mV in the female and -53.87 ± 3.21 mV in the male OA model). The resting membrane potential in OA was around 17.5% more depolarized for female and 24.5% for male model compared to respective healthy resting membrane potentials. In our OA model, current plays a prominent role in regulating membrane voltage: in both males and females with OA, increased becomes the most important contributing factor to the overall depolarization of the resting membrane potential, as revealed by sensitivity analysis. Inhibiting might therefore provide a useful approach for reversing resting membrane alterations in OA and, more broadly, could be a useful target in modulating disease pathogenesis in OA.

 inhibitors are well characterized for e.g. cardiac disease treatment, and new selective inhibitor drugs are currently in the pipeline [28]. To test whether block could be useful in overcoming OA-associated changes in resting membrane potential, we simulated block. With varying degrees of inhibition, the resting membrane potential in disease can be restored to healthy, pre-OA levels in both males and females. Therefore, re-purposing blockers feasibly provide a useful new therapeutic approach for OA treatment. blockers have, in fact, already been applied for OA treatment: clinical studies have reported positive outcomes in patients, revealing that blockers, Ouabain and Digoxin, induced a reduction in pain and elicited an apparent OA-protective efect [28], supporting our simulation results implying blockers as a potential treatment for OA.

Given this putative critical role of Na+ concentration in healthy and OA chondrocytes, novel future experiments should consider thoroughly probing Na+ balance in human chondrocytes. For instance, Figure 3.3 shows that in male and female OA chondrocytes, intracellular sodium quickly becomes depleted and remains very small, close to zero, for the rest of the simulation, suggesting that alternate sodium regulation not currently accounted for in base models may apply in OA. Model development should additionally move towards incorporating new channel species to regulate Na+ concentrations in control and in OA, including epithelial sodium channels (ENaC), which are not currently included in the model [31, 34]. Other work has also focused on voltage-gated sodium channels in chondrocytes, which may also be considered.

Naturally, the lack of available biological data from new electrophysiological and other experimental studies presents a challenge for the validation of both male and female models, both in control and in OA. Further experimental studies will be instrumental in determining the ranges of parameters for the model, as well as for simulated chondrocyte behavior evaluation. In the future, further experimental data will be helpful in further developing and expanding the current chondrocyte model. While the current study, synthetic in nature, limits its direct translation to clinical application, the present work retains utility in terms of physiological exploration and hypothesis generation.

In summary, further incorporating additional species key in the chondrocyte channelome into the model and incorporating new data from emerging experiments will permit even greater insights into the broad and nuanced role of membrane dynamics in chondrocyte function and pathophysiology.

### **3.5 Conclusions**

We have further expanded a previously published chondrocyte mathematical model enabling simulations of the osteoarthritic chondrocyte channelome. Implemented changes were based on published data, and model simulation results show that the resting membrane potential in chondrocytes becomes depolarized in OA in both male and female chondrocytes, which agrees with the experimental data acquired from males. Our sex-specifc chondrocyte model reveals diferences in resting membrane potential between males and females that could potentially contribute to the higher prevalence of OA in female patients. Finally, sensitivity analyses revealed the main currents responsible for maintaining resting membrane potential in chondrocytes and potential novel therapeutic targets for OA treatment. Our OA chondrocyte electrophysiological model provides an accessible tool for subsequent studies of OA pathogenesis and drug targeting.

3 Sex-specifc OA in Chondrocytes 43

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 4 Recapitulating Functional Heterogeneity in Electrophysiologically Active Tissues**

Meye Bloothooft, Joseph G Shuttleworth, Gabriel Neiman, Ishan Goswami, and Andrew G Edwards

**Abstract** Inter-cellular heterogeneity is central to the dynamic range and robustness of function in many tissues, particularly electrically excitable tissues. In pancreatic islet -cells, inter-cellular heterogeneity underlies the range of insulin response to glucose. In human-induced-pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs), inter-cellular heterogeneity presents a key challenge for drug screening applications. In this study, we assess the ability to reconstruct inter-cellular heterogeneity *in silico* by applying a "population of models" (PoMs) framework, i.e. collections of computational cells created via Monte Carlo variation of model parameters. We defne parameter variation based on experimentally observed heterogeneity in properties such as ion current conductances and enzymatic afnities. We then assess the accuracy of those reconstructions, based on the degree to which variation in PoM outputs (e.g. action potential duration) matches experimentally observed variation. We report that this "ground-up" approach underestimates functional heterogeneity in the hiPSC-CM population, but overestimates it in adult human cardiomyocytes. In contrast, the -cell PoM captures three distinct and physiologically relevant subclasses of -cell function. In the future, we expect PoM approaches like these will

Meye Bloothooft

Department of Medical Physiology, Division of Heart and Lungs, University Medical Centre Utrecht, NL

Joseph G Shuttleworth Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, Nottingham, UK

Gabriel Neiman Department of Bioengineering, University of California, Berkeley, USA

Ishan Goswami

Department of Bioengineering, University of California, Berkeley, USA, Department of Materials Science and Engineering, University of California, Berkeley, USA

Andrew G Edwards (corresponding author)

Department of Computational Physiology, Simula Research Laboratory, NO e-mail: andy@ simula.no

permit incorporation of realistic cellular heterogeneity in detailed models of intact tissues, and thereby aid development of sophisticated tissue-engineered platforms for therapeutics.

### **4.1 Introduction**

Functional heterogeneity among cells within a tissue is a characteristic and relevant feature of living systems. Even within a genetically homogeneous population, individual cells often show a considerable amount of phenotypic heterogeneity [1], and this feature has implications in the plasticity and robustness of the organisms. Plasticity allows adaptive cellular responses to perturbations in the environment while maintaining a robustness in functional behaviour of the tissue even after a perturbation or insult. For example, when a region of the heart becomes infarcted (dies), non-infarcted healthy areas will adapt and increase their force of contraction to maintain blood fow to the body. This ability to adapt is a form of robustness of organ function, but must be achieved in a manner that maintains the heart's electrophysiologic stability. Similarly, in obesity induced diabetes, hyperlipidemia and hyperglycemia exert excessive insulin demand on pancreatic islet -cells. In some -cells this prolonged exertion leads to exhaustion of insulin reserve and cell death. Thus, increased demand for insulin must be achieved via an adaptive response in the remaining viable -cells, which are phenotypically distinct. In both of these examples, cellular heterogeneity is central to the plasticity and robustness that allows both organs to maintain function after pathophysiological challenge. Functional heterogeneity may also broaden the range of responses to drugs [2]. For example, glibenclamide has diferential efects on pancreatic islets owing to functional heterogeneity of their -cell populations [3]. Thus, the study of cellular heterogeneity for both basic understanding of living systems and discovery of therapeutics has motivated the development of experimental techniques and the construction of computational models to simulate the dynamics and efects of cellular heterogeneity [4].

A range of single-cell experimental techniques are available to investigate this heterogeneity. These range from transcriptome assessment via single cell RNA sequencing, to classical cell electrophysiology via patch clamp. While these approaches are generally both time and labor intensive, data collected by those methods can be stored and interrogated via efcient modelling frameworks, such as the population of models approach (PoM) [5, 6]. These approaches simulate the dynamics and functional implications of heterogeneity in cell populations, including responses to pathologic insults and drug challenges.

In the PoM approach, heterogeneity is implemented by varying parameters (such as the ionic conductance of a single or multiple ion channel/s) of a baseline model, which itself captures the average characteristics of a particular cell type. Each unique confguration of parameters creates a new instance of the cell model, which is analogous to a new "member" cell of the population. In general, this procedure has been performed by frst specifying equal distributions over which all parameters of interest (e.g. all ionic conductances) can vary, and then selecting only members of that population that exhibit behaviour within the experimentally observable range [5, 6]. For example, a cardiac myocyte population could be constructed by frst allowing all ionic conductance parameters to be sampled from Gaussian distributions with an arbitrary standard deviation of 30% of each mean parameter value in the baseline model - this sampling strategy involves homogeneous variability across parameters. Some members of this PoM will exhibit parameter combinations that are biologically unrealistic, and thus possibly exhibit behaviour outside the experimentally observable range. To eliminate those models, it is common to reject models that produce outputs or features, such as a cardiomyocyte action potential duration, that are outside the range of experimental observations for the respective features. Continuing the example, a cardiomyocyte population member will be deemed non-viable if the action potential duration is far longer or shorter than observed experimentally. Similarly, a non-fring -cell model will be rejected from the islet PoM. Because, in these examples, the PoMs are generated from arbitrary and homogeneous parameter distributions, they will exhibit appropriate behavioural characteristics, but may not involve truly biologically representative parameter combinations. Thus, the behaviours may not result from truly representative underlying dynamics.

In this work, we apply a diferent approach, with the intention of introducing more biologically realistic variability in PoMs of three diferent cell types: human adult cardiomyocytes, hiPSC-CMs and primary islet -cells. By introducing variability to the input parameters that is grounded in *experimentally observed variability for each parameter*, we hope to span the biologically reachable confgurations and limit the number of unrealistic confgurations that are generated. That is, we use published experimental variation (standard deviations) for ionic conductance parameters and enzymatic activity/afnity parameters in order to apply parameter-specifc heterogeneity rather than variation that is homogeneous across the varied parameters. Based on these constraining conditions, we demonstrate the feasibility of PoMs to generate a heterogeneous population representing diferent functional phenotypes.

### **4.2 Methods**

An illustrative example of our approach in generating PoMs is shown in Figure 4.1. To recapitulate the functional heterogeneity in the tissue PoMs, we introduced parameter-specifc variability into our models. To do this, we identifed a set of parameters which we expect to account for much of the inherent biological variability infuencing electrophysiologic function in these cell types. These parameters are a selection of maximal conductance (or fux) parameters for both the cardiomyocyte models (hiPSC-CM and adult) and the -cell model. Additionally for the -cell, a selection of kinetic parameters for enzyme reactions in the glycolysis cycle were varied to introduce metabolic heterogeneity. We chose to endow each of these pa-

Fig. 4.1: Illustrative example of generating population of model (PoM) of human cardiomyocytes and islet -cells. Experimental data obtained from literature was used to generate log-normal distributions of electrophysiological and/or metabolic parameters used as inputs to computational model. Sets of features are obtained from the simulations that are then used to generate viable candidates for the PoMs representative of cardiomyocyte or -cells via constraining conditions.

rameters, , with a log-normal distribution, such that log ∼ ¯ (0,<sup>2</sup> ) where ¯ is the original, unmodifed parameter from the baseline model and is a standard deviation derived from the literature. We term the random variable ¯ the *scaling factor*. To assign values for the various , we surveyed the literature for measured variability of each parameter of interest. To mitigate the impact of systematic diferences among the means of diferent datasets, and to keep our parameter distributions centred on the original parameter values, we computed the coefcient of variation for each parameter ():

$$
\sigma\_p \coloneqq \frac{\sigma\_{\text{lit}}}{\mu\_{\text{lit}}} \tag{4.1}
$$

where lit and lit are the mean and standard deviation, respectively, of the observations found in the literature.

The published mathematical models chosen as the baselines for the adult-CM and the hiPSC-CM were developed by [7] and [8], respectively. For the -cell PoM, the electrophysiological model developed in [9] and glycolysis model developed in [10] were used to generate the heterogeneous population of -cells. In the following subsections, we will detail the generation of parametric distributions and the constraining conditions via feature extraction for each PoM.

### **4.2.1 Cardiomyocytes**

The O'Hara-Rudy model (ORd, [7]) is a widely-used mathematical description of the average human left ventricular cardiomyocyte. We specifcally used the epicardial version for our adult human baseline model. Likewise, the baseline model for hiPSC-CMs was that published in [8]. Both models use the following fundamental equation for calculation of time-varying membrane potential:

$$C\_m \frac{dV}{dt} = -\sum\_i I\_i \tag{4.2}$$

where is the capacitance of the cell membrane and each is a current describing the transmembrane current carried by each species of ion channel present in the model. For each ion channel, there is a maximal conductance parameter , or analogous permeability parameter, such that = ˆ for some ˆ . ˆ is in turn described by a system of ordinary diferential equations that specify the voltage, ion concentration and time-dependencies of ion channel function. For all currents, we only vary the maximal conductance (), and not the parameters of the remaining equations contributing to ˆ . The formulation for sarcoplasmic reticulum (SR) calcium reuptake is diferent to those of the ion channels, but an analogous term that refects the maximal rate of SR calcium reuptake ¯ SERCA. As for the ion channel conductances, this parameter was varied to introduce heterogeneity in this key element of intracellular calcium handling.


Table 4.1: **Scaling parameters for the cardiomyocyte PoM**

Sources to construct coefcient of variation () listed by the reference numbers

[11], <sup>2</sup> [12], <sup>3</sup> [13], <sup>4</sup> [14], <sup>5</sup> [15], <sup>6</sup> [16] [17], <sup>8</sup> [18], <sup>9</sup> [19], <sup>10</sup>[20], <sup>11</sup>[21], <sup>12</sup>[22], <sup>13</sup>[23], <sup>14</sup>[24], <sup>15</sup>[25], <sup>16</sup>[26], <sup>17</sup>[27], <sup>18</sup>[28], <sup>19</sup>[29], <sup>20</sup>[30], <sup>21</sup>[31], <sup>22</sup>[32]

Fig. 4.2: A. The resultant log-normal distributions for the Ks scaling factor, along with histograms resulting from 100 samples, are shown for the representative parametric variation in adult CM and the hiPSC-CM inputs for the models. B. The scaling factors that were sampled in order to create the PoM are shown for two cardiomyocyte models. Parameters from the models that were excluded from the POM are shown in grey, whereas the parameters which were calibrated into the PoMs are highlighted in orange.

Table 4.1 provides the values for each varied parameter for both the adult human myocyte model and the hiPSC-CM model. The variability presented in this table may not completely capture the variability in each maximal conductance parameter. Nevertheless, these values provide a reasonable approximation, and yield plausible distributions for each of the relevant parameters. We took 500 samples from each of these distributions and the results are shown in Figure 4.2A alongside the relevant sampling distributions.

The PoMs generated by varying the input parameters as summarized in Table 4.1 were solved via MATLAB stif-solver ode15s, and allowed to equilibrate by simulating for 500s before computing the output features. Calibrated PoMs for both the adult CM and hiPSC-CM were obtained by rejecting models that produced features outside the experimentally observed range. We identifed a collection of features from the literature: *maximum upstroke velocity (voltage)*, *minimum diastolic potential*, *voltage amplitude*, *action potential duration (APD*90*)*, *maximal departure velocity for Ca2+*, *diastolic Ca2+ concentration*, *Ca2+ transient amplitude*, *Ca2+ transient time-to-peak*, and *time constant of Ca2+ decay*. Representative scaling factors sampled from the output of the models obtained via the constraining criterion are shown in Figure 4.2B. In the results section, we discuss the typical features obtained and constraining criteria to generate calibrated PoMs for the cardiomyocytes.

### **4.2.2 Islet -cells**

A robust secretion of insulin by the -cell upon glucose challenge relies heavily on the coupling of the metabolic oscillations in the glycolysis pathway with electrical oscillations (e.g. isolated action potentials and action potential bursting). Recognizing this, the -cell PoM was generated via a coupling of the electrophysiological model developed in [9] and glycolysis model developed in [10]. A coupling of these models was demonstrated in [33], albeit using single mean values for the electrical and glycolytic component parameters. A simplistic illustration of this coupled model is shown in Figure 4.3. To delineate contributions made by the glycolytic and electrophysiologic systems to the overall functional heterogeneity in our -cell populations, we divide the parametric distributions into two components: ionic and glycolytic.

Fig. 4.3: Schematic representation of a -cell's metabolic and electrophysiological circuit involved during the glucose-stimulated-insulin-secretion. Included in the schematic are representative formulae of the glycolytic component's enzyme rate equations and a few parameters (in red) that are changed to generate PoMs. Lognormal distributions for two of the 14 parameters are shown in the bottom right hand side of the fgure.

Details of the ionic models can be found in the cited articles, but briefy electrical oscillations in the system were modeled as described by Equation 4.2. The maximal conductance () of 8 ion channels were varied based on experimental data in the literature (references found in [34], [35], and [36]). These variations, reported via coefcients of variation () and their sources, are summarized in Table 4.2.

The glycolysis model comprised enzymatic reactions from glucokinase to glyceraldehyde-3-phosphate. The rate of ATP generation through these steps and via downstream oxidative phosphorylation in the mitochondria is introduced in a phenomenological manner via the variable, *a*, which represents ATP level.

$$\frac{da}{dt} = V\_{GAPDH} - k\_A a \tag{4.3}$$

where is the metabolic fux due to GAPDH and is a phenomenological time constant.The ion channel conductance of channel is then altered via:

$$\text{g}\_{KATP} = \frac{\text{\hat{g}}\_{KATP}}{1 + a} \tag{4.4}$$

In our study, we varied the kinetics of the enzymatic reactions by varying the limiting rates , of glucokinase, phosphofructokinase, and glyceraldehyde 3-P dehydrogenase. In addition we varied other kinetic parameters such as the halfactivation concentrations 0.5 . These parameters are highlighted in red in our Figure 4.3. Variations in these glycolytic components of our model and their sources are summarized in Table 4.2. This variability was again assumed to be log-normal distributed and implemented as such for construction of our PoMs. Sources for these datasets can be found in [37], [38],[39],[40].

Since the behaviour of individual -cells is much more diverse than those of cardiac myocytes, the constraining conditions are also less well-established. Thus, we discuss the features used for constraining the -cell PoMs in the results, and we report sub-classes of -cells obtained via our approach. Furthermore, we will analyze the efects of perturbations to the glycolysis and ionic components in maintaining functional heterogeneity.

### **4.3 Results**

### **4.3.1 Cardiomyocyte PoM**

As described in the methods, we allowed each PoM member model to equilibrate for 500s before computing the output features. The last 10 action potentials of each simulation for these models were used to extract features that were compared against the constraining criteria to generate the calibrated PoMs. Constraining criteria were:

#### 4 Functional Heterogeneity 53

### Table 4.2: **Scaling parameters for islet -cell PoM**

Sources to construct coefcient of variation () listed by the reference numbers


[34], <sup>2</sup> [35], <sup>3</sup> [36], <sup>4</sup> [37], <sup>5</sup> [38], <sup>6</sup> [39], <sup>7</sup> [40]


The full populations, showing both the accepted (orange) and discarded (grey) model confgurations are shown in Figure 4.4A, together with the percent of confgurations discarded based on each pairwise combination of constraints in Figure 4.4B. This fgure provides a frst indication of how well variability in the model behaviour is captured by applying a data-defned variation in the model inputs. A large number of confgurations had to be discarded for both models, although the reasons for exclusion were diferent. Specifcally, we retained only 55 of 500 models for our adult ventricular cardiomyocyte PoM, and 60 of 100 models for our hiPSC-CM PoM. For the adult ventricular PoM the vast majority of exclusions were made because the confguration failed either <sup>1</sup> (43% of confgurations) or <sup>9</sup> (40% of confgurations). This can also be seen in Figure 4.4A, where many confgurations fail to repolarize below -65 mV (9), whereas others reach supraphysiologic membrane potentials (either > 70mV or < −100mV). For the hiPSC-CM model the reasons for exclusion were physiologically reachable but not observed in published hiPSC-CM phenotypes. That is, the vast majority of excluded confgurations simultaneously failed criteria <sup>2</sup> - 9, each of which was chosen to address a diferent source of biologically reachable instability. This can also be seen in Figure 4.4A, where most of the discarded (grey models) have attracted to a second equilibrium in membrane potential between −10mV and −20mV. This property of electrophysiologic bistability is a real biological property of cardiac myocytes [41]. It is thus a fair question to ask whether such confgurations should be discarded for the hiPSC-CM PoM, as they may in fact be realistic members of most hiPSC-CM cell populations. Regardless, this initial characterization suggests that the experimentally-defned variation in hiPSC-CM input parameters results in sharp transition to gross electrophysiologic instabilities, most of which remain physiologically plausible. On the other hand, the majority of adult ventricular PoM exclusions resulted from either biologically unreachable phenotypes (V > 70mV or < −100mV), or to an unstable resting potential (minimum V > −65mV).

Fig. 4.4: **A** The action potential and calcium transient phenotypes from the adult CM and hiPSC-CM PoMs. Traces of cells accepted into the calibrated PoMs are in orange, while those discarded are shown in grey. **B** Tables summarizing the % of models excluded by the constraining criteria. Diagonal values show the percentage of models excluded by each criterion in isolation, whereas the of-diagonal entries show the % of confgurations for which each combination of two excluding criteria were present. The meaning of each criterion, i.e. , is explained in the main text.

Figure 4.5 shows the relationships between major features of the fnal calibrated PoMs for the two cardiac models, and provides a basis for comparing systematic diferences in the phenotypic span of each PoM. The stronger relationships for the hiPSC-CM PoM between Ca2<sup>+</sup> handling features and AP amplitude indicate the larger role of sarcolemmal fuxes in determining Ca2<sup>+</sup> features. In contrast, the adult

Fig. 4.5: Scatter plots of feature pairs from both PoMs. The ftted distributions are shown along the diagonals using a kernel density estimator. Noteworthy diferences between the two PoMs are boxed in the red (adult) and purple (hiPSC-CM).

ventricular PoM Ca2<sup>+</sup> dynamics are more internally determined, as demonstrated by stronger relationships boxed in green. These diferences refect intrinsic characteristics of the baseline hiPSC-CM and adult ventricular cardiac myocyte models. They are a well-known diference in the biological underpinnings of the two cell types and refect the generally less "mature" Ca2<sup>+</sup> handling phenotype of hiPSC-CMs relative to adult ventricular myocytes.

Finally for the cardiac PoMs, Figure 4.6 shows important diferences between the two calibrated PoMs in terms of degree to which the variability in their AP and Ca2<sup>+</sup> handling features refect those reported in the literature. Specifcally, it is clear that the adult ventricular PoM features (top row, orange histograms) are generally more variable than is reported experimentally (dashed lines). Only the APD<sup>90</sup> feature for this PoM appears to approximate reported experimental variation. In contrast, the hiPSC-CM feature distributions (bottom row, orange histograms) are less variable than published reports of these phenotypes (dashed lines). This presents a fundamental distinction, and suggests that literature reports of variation in the underlying currents and fuxes overestimate the true physiologic variability in human adult cardiac ventricular myocytes. In contrast, the same techniques applied to hiPSC-CMs underestimate observable phenotypic variability. We explore this important observation in the discussion.

CaT amplitude / umolL <sup>1</sup> APD90% / s Fig. 4.6: Histograms showing the variability in the features of the PoM. The frst row shows the adult ventricular cardiomyocyte PoM, and the second row shows the hiPSC-CM PoM. In each plot, the probability density function of a representative Gaussian distribution shows the variation typically found in measurements of these features in the literature. These distributions have been normalised such that they are centered on the means of each histogram.

### **4.3.2 Islet -cell PoM**

The individual models for the -cell PoM were simulated to capture cell behaviour over 10 minutes - the acute phase of glucose-stimulated insulin release. The parameter variations provided in Table 4.2 were used to generate 1000 cells to assess the importance of 3 sources of heterogeneity, by varying model parameters associated with: (a) only glycolytic (enzyme) model components, (b) only ionic (ion channel and transporter) components, (c) both ionic and glycolytic components. In all three cases, each model confguration was simulated under high glucose (10 mM) and low glucose (2 mM) conditions. Thus, the total number of simulations across all 3 cases (3000 cell model confgurations), at both glucose concentrations, was 6000.

Figure 4.7A shows V and cytosolic calcium for representative simulations at high glucose. We see cells 1 through 3 demonstrate similar features as experimentally observable sub-classes of -cells, i.e. spiking, bursting, and plateau phenotypes. Nonfring cells such as those shown in cells 4 and 5 were also observed. As for the cardiac PoMs, in order to generate a calibrated -cell PoM, we have to invoke and apply a set of constraining criteria to reject confgurations that produce features outside the ranges of experimental observations. However, because the -cell phenotype is the most heterogeneous of all cell types assessed here, there exists far fewer such objective constraints in the -cell literature. However, some benchmarks are available. For example, [42] suggested two metrics of activity based on -cell V

Fig. 4.7: Representative simulation outputs and categorizations. **A.** Representative simulation outputs of membrane potential and cytosolic calcium representing different sub-classes of -cells. 1. Spiking cells 2. Bursting cells 3. Plateau cells 4 & 5. Non-fring cells. Scales on the plots represents time of 240s. **B.** Plot of activity fraction vs. Δpeak for the 5 representative cells. **C.** Number of peaks in the Fast Fourier Transform (FFT) of the membrane potential traces for the 5 representative cells.

(the activity fraction and Δpeak), as means of classifying phenotypic sub-classes. Activity fraction was defned as the fraction of time V was above an arbitrary threshold level. The defnition of Δpeak was set as the diference of the two major components of V-time probability distribution (i.e. the diference between "resting" V and "active" V). As for [42], we did not observe clear delineation of the three sub-classes, nor could we readily diferentiate fring vs. non-fring -cells via these metrics (Figure 4.7B). However, we surmised that Fast Fourier Transform (FFT) of the V signal may allow more clear distinction of these classes. Figure 4.7C shows the number of FFT peaks corresponded with the number of dominant harmonics in the V signal for each representative recording in Figure 4.7A.

Based on this analysis, we set the constraining criteria as follows. Any cell with less than 100 peaks in the V peaks, was considered to be quiescent and thus discarded. For cells that had more than 100 peaks in their voltage signal, FFT was performed on the voltage signal to classify them into spiking, bursting, or plateau cells. Based on our initial analyses, cells with less than 20,000 peaks in the voltage FFT were classifed as plateau cells. Bursting cells were classifed as those having anywhere above 20,000 FFT peaks but less than 40,000 peaks. Finally, cells with FFT peaks greater than 40,000 were classifed as spiking cells.

Fig. 4.8: **A.** Percentage of viable cells based on the chosen constraining conditions for six diferent cases (1000 runs each). **B.** Distribution of the sub-classes of -cell from the viable candidates for each of the six cases. HG: High glucose (10 mM). LG: Low glucose (2 mM). Analyses of the cytosolic calcium traces of each of the sub-classes of -cells in the viable population obtained with high glucose: **C.** mean values, **D.** peak values, **E.** Number of FFT peaks. Statistical diferences between the groups were determined by one-way ANOVA followed by Tukey's HSD post-test.

Based on these constraining criteria, we quantifed the viable cells that were not rejected for both glucose concentrations and all 3 cases of heterogeneity (glycolytic, ionic, both), as shown in Figure 4.8A. The greatest reduction in viable cells was observed for glycolytic heterogeneity only, thus highlighting the narrow metabolic parameter range functional -cells operate under and its importance for achieving robust insulin secretion. Another interesting aspect of these data was the number of active cells at low ("resting") glucose (2 mM). This feature of basal activity is unique to human -cells, and underlies basal insulin release at the lower blood glucose concentrations present in humans. Thus, our PoM approach was able to capture this aspect of human-specifc -cell physiology.

Further analysis of the viable candidates allowed us to assess the fraction of each of the 3 PoMs belonging to the sub-classes of -cell phenotypes at each glucose concentration (Figure 4.8B). The lowest fraction of bursting cells was observed in the low-glucose (LG) cases when both components or only ionic components were varied. Upon altering glycolytic components alone, plateau cells comprised the only viable sub-class, thus suggesting the importance of coupling the ionic and glycolytic elements to generate the full range of functional phenotypes observable in human -cell populations. Mean cytosolic calcium of the viable cells within our PoM were not statistically diferent between spiking and bursting cells, suggestive of less pronounced basal shift of calcium signalling (Figure 4.8C). However, plateau cells exhibited a mean cytosolic calcium signature distinct from the other two classes, while peak cytosolic calcium grouped plateau and bursting phenotypes. (Figure 4.8D). Applying FFT analysis to the calcium signals also showed that the number of harmonics in the signal could distinguish the sub-classes of -cells, thus suggesting that this may be a viable means of distinguishing functional heterogeneity in experiments, for which calcium recordings are considerably easier.

### **4.4 Discussion**

In this manuscript, we report the construction of PoMs for three electrophysiologically active tissues that were calibrated against experimental data from the literature. PoMs have proven useful in investigation of cardiac electrophysiological variability and recent studies have furthered the methodology by explicitly incorporating experimental data into the construction of populations of models [43, 6]. In our CM PoMs, we observed that the fuctuations in the input parameters for the models allowed us to generate cells with features comparable to those observed in both primary and hiPSC-derived CMs, albeit with difering coverage of the range of function observable in experiments. In both cardiac PoMs we assumed the main source of variability is difering expression of ion transporters, and that this variation is responsible for the observable variation in their V and calcium cycling features. In hiPSC-CMs particularly, this variation is likely to be pronounced given their immature ("fetal-like") phenotype and the degree to which this can be matured. Importantly, we observed that experimentally-defned variation in the input parameters of the cardiac models resulted in functional output that difered in its range for the two models. That is, for the adult CM PoM, the range of AP and calcium handling features was broader than observed experimentally, whereas it was narrower than experimental for the hiPSC-CM PoM. This is a fundamental observation about the agreement between two types of data recorded in these cell types. For the adult CM it could be that biological covariance (which we have not applied in PoM construction) among the ionic currents reduces the range of reachable AP and calcium features in real cells, or that the relative scarcity of adult human CM data means that these two classes of data are not yet fully internally consistent. In contrast, the greater experimental variability for hiPSC-CM phenotype features suggests either that experimental ionic current measurements are subjected to overly constraining inclusion criteria. Alternatively, our hiPSC-CM model formalism may somehow be constrained to be overly stable with respect to the phenotypic variation resulting from variation in the input currents. Regardless, these questions are fundamental to understanding how perturbations to either adult or hiPSC CMs should be expected to impact those tissues *in vivo*.

The islet -cell PoMs were generated via alterations in both the parameters involved in glycolysis model components as well as the ionic conductance components. We observed that the number of harmonics in V of the simulated computational cells was enough to distinguish diferent sub-classes of -cells observed experimentally, although we believe the proposed constraining criteria requires further validation via prospective screening of experimentally obtained data sets. Furthermore, we observed that the functional heterogeneity of the -cells is dependent on the coupling of the metabolic and ionic components, as noted by the number of viable cells and their sub-types when only glycolytic components were varied. This, in part, explains why glycolytic bottlenecks in stem cell-derived -cells show non-robust glucose-stimulated-insulin-secretion behaviour [44, 45]. The calibrated PoMs generated in this study were able to capture electrical oscillations at low glucose, which is a unique feature of human -cells (compared to rodents). Thus, the methodology developed should be applicable for creating human-specifc 3D models of coupled -cell clusters with and without partner cells such as -, -, - and pancreatic polypeptide cells. When coupled with experimental studies, these 3D models could be invaluable to understand the efect of heterogeneity on intact islet function, and to develop sophisticated, robust tissue-engineered platforms for therapeutics.

Here we have demonstrated the use of experimental constraints to frst construct and then calibrate PoMs of three human cell types. The cardiomyocyte PoMs partially reproduced observable variation in the functional features for their respective cell types, albeit with clear systematic inaccuracies. The ability of the human -cell PoM to recapitulate the range of function observed in experimental -cell recordings is encouraging for applying this strategy to integrated function of -cells in human islets. Constructing new methods for constraining these PoM approaches will require continued advancement of single cell assays.

4 Functional Heterogeneity 61

### **References**


pluripotent stem cell derived cardiomyocytes: a dynamic clamp study with virtual IK1. *Frontiers in Physiology*, 6, 2015.


stem cell-derived β cells in vitro is inhibited by a bottleneck in glycolysis. *Cell Reports*, 31(6):107623, May 2020.

45. Ishan Goswami, Eleonora de Klerk, Phichitpol Carnese, Matthias Hebrok, and Kevin E Healy. Multiplexed microfuidic platform for stem-cell derived pancreatic islet cells. *Lab on a Chip*, 2022.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 5 Realizing Synaptic Signal Transmission During Astrocyte-Neuron Interactions within the EMI Framework**

Julia Gorman, Konstantin Holzhausen, Joyce Reimer, and Jørgen Riseth

**Abstract** The tripartite synapse or "neural threesome" refers to the interplay in the synapse between neighbouring neurons, the synaptic cleft, and the surrounding glial cells. Despite extensive research, the efects of glial cells, such as astrocytes, on signal transduction between neurons are not fully understood. The Kirchhof-Nernst-Planck (KNP) and Extracellular-Membrane-Intracellular (EMI) models constitute a promising framework for modeling these kinds of systems. However, they lack the neurotransmitter-related mechanisms that are necessary to bridge signal transduction across the synaptic cleft. Here, we propose an extension to the KNP-EMI model by a spatio-temporal difusion-based description of the most prominent neurotransmitter, glutamate, that allows for investigation of the contribution of astrocytes to the functionality of the synapse. We validate our model by showing that the presence of an astrocyte in the domain afects the glutamate fux across the postsynaptic terminal, as observed physiologically. The proposed extension ofers a sufciently simple way of integrating synaptic glutamate dynamics into the KNP-EMI framework. It introduces the relevant interactions between electrical activity and difusion processes at the tripartite synapse that are necessary to assess how astrocytes might contribute to the functionality of the synapse. This work has implications for future studies involving glial mechanisms and other charged species within the KNP-EMI framework.

Jørgen Riseth (corresponding author)

Julia Gorman Department of Neurosciences, University of California, San Diego, USA

Konstantin Holzhausen Department of Physics, University of Oslo, NO

Joyce Reimer Division of Biomedical Engineering, University of Saskatchewan, CA

Department of Numerical Analysis and Scientifc Computing, Simula Research Laboratory, NO e-mail: jorgennr@simula.no

### **5.1 Introduction**

Circuits in the brain are composed of complex connections of both neurons and glia. Neurons communicate via patterns of electrical signals that are propagated by molecules and ions. These signals, called action potentials, communicate information that gives rise to cognitive function, behaviors, and movements. Action potentials travel from neuron to neuron across small gaps called synapses. The neuron sending the signal, the presynaptic neuron, releases specifc molecules called neurotransmitters to the receiving neuron, or the postsynaptic neuron. However, neurons are not the only cells that are found in the synapse. A type of glial cell called an astrocyte envelopes the synapse and alters information transmission by modulating the neurotransmitters as they difuse across the synapse [1]. Thus, together, the neurons, the synapse, and the astrocyte form a "neural threesome" [2], or tripartite synapse. The role of this astrocyte mediated modulation in electrical and chemical interplay is not well understood.

One of the ways astrocytes are thought to afect synaptic transmission is through reuptake of the neurotransmitter glutamate. Once an electrical signal reaches the end of the presynaptic neuron, prepackaged vesicles release glutamate into the synaptic cleft. After release, glutamate difuses across the synapse where it binds AMPA receptors on the postsynaptic neuron, resulting in a net inward sodium (Na<sup>+</sup> ) current through the ligand-gated AMPA receptor ion channels. As a result, the intracellular potential starts to increase which triggers the opening of nearby voltage-gated Na<sup>+</sup> channels. This channel opening allows for a larger Na<sup>+</sup> infux, bringing the postsynaptic neuron closer to its threshold for depolarization. If the signal is strong enough, subsequent ion dynamics across the postsynaptic membranes result in the onset of a new action potential and signal propagation continues.

As glutamate difuses across the synapse, astrocytes take up a percentage of the glutamate ions through high-afnity glutamate transporters [3, 4]. This mechanism is suggested to be a way of reducing cellular toxicity, as too much glutamate in the synapse can have excitotoxic efects [5]. Recently, it has been found that astrocytes also release glutamate to neurons in response to extracellular glutamate cues [6]. Because of this dual mechanism of glutamate uptake and release, astrocytes act as a bufer to balance excitatory and inhibitory neuronal activity.

Here, we propose a simple mechanistic model of glutamate in the tripartite synapse that allows for the assessment of these questions within the Kirchhof-Nernst-Planck (KNP) - Extracellular-Membrane-Intracellular (EMI) modeling framework. The KNP-EMI model is a mechanistic mathematical model for excitable tissue in spatiotemporal resolution. It is stated in the form of a set of coupled partial diferential equations, and includes interactions between ionic transmembrane currents and electric potentials [7, 8, 9]. In [7], Ellingsrud *et al.* show that their KNP-EMI model accurately describes depolarization along excitable membranes. Thus, it constitutes a truly spatio-temporal model of signal transduction along a cell, for example, along an axon or a dendrite of a neuron. However, the model in its current form does not account for neurotransmitters in the extracellular space, nor for neurotransmitterrelated mechanisms and interactions. For this reason, the KNP-EMI model is not able to capture the regulation of signal transduction at the synapse.

We propose an extension by introducing glutamate to the model. In order to bridge the synaptic gap, we identify three relevant biophysical mechanisms: glutamate release into the cleft by exocytosis from the presynaptic neuron, difusion in the extracellular space, and binding to postsynaptic AMPA receptors. All three of these processes are considered in our implementation. In addition to this, we assume that glutamate uptake is the fundamental interaction of the astrocyte at the synapse whereby astrocytes regulate signal transmission. Therefore, our proposed glutamate difusion model describes the properties of signal transduction unique to the synapse; in particular, how astrocytes afect the glutamate gradient and how a changing glutamate gradient afects signal propagation. Our extension to the KNP-EMI model features a simplifed yet qualitatively complete description of signal transmission across neurons. The extended model will allow for a comprehensive investigation of the astrocyte's role in signal propagation across neurons in a way that takes spatial efects, such as explicit geometries and local concentrations, into account.

The EMI model constitutes a general framework for modeling excitable tissue ranging from cardiac cells to neural systems [10]. Besides its broad range of application, it allows for the incorporation of explicit morphologies within the system of interest. In combination with the KNP framework, the KNP-EMI model allows for the representation of a complicated geometry and detailed electrochemical profles [8]. In order to leverage its power, we design our glutamate model under consideration of its integrability into the KNP-EMI model. Although its formulation is independent of the respective morphology, our implementation assumes a simplistic geometry in order to show its functionality as a proof of concept.

Our paper is structured as follows. First, we familiarize the reader with the geometry of the problem conceptually and derive the glutamate model. Subsequently, we demonstrate a way to estimate the model parameters based on literature values, restricting our model to hippocampal synapses. After providing implementation details, we validate the glutamate model by showing that it behaves as intended. Finally, we discuss the implications of our model within the context of the KNP-EMI framework.

### **5.2 Methods**

The glutamate model consists of two parts. The frst part describes the transport of glutamate across the synaptic cleft by difusion, whereas the second part describes how AMPA receptors react in response to glutamate to trigger a Na<sup>+</sup> current across the postsynaptic terminal. Here, we describe the computational domain and the governing equations for the glutamate model. We then describe the parameters that are used within the model, and fnally, the numerical methods used in our simulations.

### **5.2.1 Representation of the Computational Domain**

The tripartite synapse features a variety of membranes that behave diferently. The characteristic behavior of each membrane is determined by the density and type of functional components it exhibits, e.g. channels, pumps, and receptors. Figure 5.1 shows a schematic of our envisioned computational domain for the full KNP-EMI model. It resembles a simplifed and abstracted two-dimensional realization of a tripartite synapse setup. The domain is divided into two distinct subdomains: the

Fig. 5.1: Representation of the computational domain in our glutamate KNP-EMI model. Synthesis of computational domains, channel dynamics, and boundary conditions.

extracellular space, Ωe, and the intracellular space, Ω<sup>i</sup> . Ω<sup>i</sup> represents the interior of the postsynaptic neuron's dendrite. The interface boundary separating Ω<sup>e</sup> from Ω<sup>i</sup> consists of the postsynaptic terminal, Γpost, and extrasynaptic membrane of the dendrite, Γd. We assume that Γ<sup>d</sup> is impermeable to glutamate. In contrast, Γpost is a region with an afnity for glutamate due to the presence of AMPA receptors, which are only located on the part of the postsynaptic membrane that faces the synapse [11]. The exterior boundary consists of three biophysically distinct types. Γpre denotes the membrane between the presynaptic neuron and the synapse. Γ<sup>a</sup> represents the membrane on the astrocyte, which also has an afnity for glutamate due to the presence of glutamate transporters that clear it from the synapse [4]. The remaining part of the exterior boundary is open for glutamate to pass, since the neurotransmitter freely difuses through the extracellular space out of the synaptic region. We note that because this report only covers the glutamate model, our computational domain is for now restricted to the extracellular space, Ω.

### **5.2.2 Mathematical Modeling of Glutamate Dynamics**

The glutamate transport is modelled using a pure difusion equation: fnd the glutamate concentration = (,) such that

$$
\partial\_t \lg \mathbf{g} = D\_\mathbf{e} \Lambda \mathbf{g} \quad \text{on} \quad \mathfrak{Q}\_\mathbf{e} \times [0, T] \quad . \tag{5.1}
$$

Here, is the difusion coefcient for glutamate within the extracellular space. The glutamate release from the presynaptic terminal and absorption of glutamate on the various cell membranes are modelled using the following combination of von Neumann and Robin boundary conditions:

$$-D\_{\mathbf{e}} \nabla \mathbf{g} \cdot \mathbf{n}\_{\Gamma\_{\mathbf{d}}} = \qquad \qquad \mathbf{0} \qquad \qquad \text{on} \quad \Gamma\_{\mathbf{d}} \times [0, T] \tag{5.2}$$

$$-D\_{\mathbf{e}} \nabla \mathbf{g} \cdot \mathbf{n}\_{\Gamma\_{\mathbf{a}}} = \begin{array}{c} k\_{\text{update}} \ \mathbf{g} \vert\_{\Gamma\_{\mathbf{a}}} \qquad \text{on} \quad \Gamma\_{\mathbf{a}} \times \begin{bmatrix} 0,T \end{bmatrix} \tag{5.3}$$

$$-D\_{\rm e} \nabla \mathbf{g} \cdot \mathbf{n}\_{\Gamma\_{\rm post}} = \quad k\_{\rm bind} \ \mathbf{g} \vert\_{\Gamma\_{\rm post}} \qquad \text{on} \quad \Gamma\_{\rm post} \times \{0, T\} \tag{5.4}$$

$$-D\_{\rm e} \nabla g \cdot \mathbf{n}\_{\Gamma\_{\rm pre}} = -\sum\_{i} v\_{\nu}^{i} \,\delta\left(t - t\_{i}\right) \quad \text{on} \quad \Gamma\_{\rm pre} \times \{0, T\}. \tag{5.5}$$

Note that here, **n** denotes the outward pointing normal vectors of the respective boundary surfaces. As a consequence, the left hand sides denote the outward fuxes of glutamate through them. We model glutamate binding to the AMPA receptors as well as glutamate uptake at the astrocyte's membrane by frst order reaction kinetics. Thus, bind specifes the afnity of the AMPA receptor towards glutamate and uptake the absorption rate. The refecting boundary condition at Γ<sup>d</sup> in Equation (5.2) realizes the impermeability of glutamate through the extrasynaptic membranes of the dendrite. Equation (5.5) describes the source of glutamate in our model, which is vesicular release due to exocytosis. In accordance with Clements *et al.* [12] and Jonas *et al.* [13], we model an instantaneous release of vesicles at release times . denotes the vesicular spatial profle defned in Equation (5.9). On the remaining open boundaries, we assume homogeneous Dirichlet boundary conditions

$$\|\mathbf{g}\|\_{\Gamma\_{\text{open}}} = \mathbf{0} \quad . \tag{5.6}$$

This choice is made for simplicity, and future work will beneft from more realistic modelling of free difusion across the open boundaries. Since glutamate is supposed to be cleared from the synaptic cleft after exocytosis, we assume an initial equilibrium of (, = 0) = 0.

### **5.2.3 Modeling AMPA Gating Dynamics**

We model the dynamics of the AMPA receptors on the postsynaptic terminal using a model proposed by Tewari and Majurmdar [14]. The ODE model describing the AMPA gating dynamics reads

$$\frac{d}{dt}\,m\_{\rm AMPA} = \alpha\_{\rm AMPA} \, g \Big|\_{\Gamma\_{\rm paa}} \left(1 - m\_{\rm AMPA}\right) - \beta\_{\rm AMPA} \, m\_{\rm AMPA} \quad . \tag{5.7}$$

Here, the glutamate concentration, , couples linearly to the opening rate of the AMPA receptor, AMPA. The ODE model governs the dynamics of the gating variable, AMPA, point-wise along the boundary Γpost, and depends on the local glutamate concentration, |Γpost . We intend to initialize the system in its equilibrium state. The initial state (, = 0) = 0 corresponds to the equilibrium state of the glutamate model. For = 0, we get ∗ AMPA = 0 as the stationary state of the gating variable ∗ AMPA at <sup>=</sup> 0. Therefore, we choose AMPA( <sup>=</sup> <sup>0</sup>) <sup>=</sup> ∗ AMPA = 0.

The gating variable, AMPA, controls an AMPA-specifc Na<sup>+</sup> current, AMPA, through the postsynaptic terminal. This current is given by

$$I\_{\rm AMPA} = \mathcal{g}\_{\rm AMPA} m\_{\rm AMPA} \left( V\_{\rm post} - V\_{\rm AMPA} \right) \quad \text{on} \quad \Gamma\_{\rm post} \tag{5.8}$$

[15, 14], with post being the local membrane potential at the postsynaptic neuron and AMPA referring to the reversal potential. In the KNP-EMI model, AMPA is associated with the Na<sup>+</sup> -specifc Nernst potential.

### **5.2.4 Estimation of Model Parameters**

All parameters for the AMPA receptors were chosen in the same way as the original model by Tewari *et al.* [14]. Following Clements *et al.* and Tewari *et al.*, we assume vesicles with a diameter Ves = 40 nm, containing glutamate at concentrations of Gl Ves = 60mM. In addition, we assume a homogeneous glutamate surface density Gl Ves across the vesicle's cross section on the presynaptic terminal. We may now defne a vesicular profle at Γpre, , as

$$v\_{\nu}^{i}(\mathbf{y}) = \begin{cases} \zeta\_{\text{Ves}}^{\text{Gl}} & \forall \, |\mathbf{y} - \mathbf{y}\_{\text{Ves}}| \le d\_{\text{Ves}} \\ \mathbf{0} & \text{else}, \end{cases} \tag{5.9}$$

where Ves denotes the center position of the vesicle on the membrane and denotes the coordinate parameterizing the one-dimensional boundary surface. We assume that the arrival of the action potential at the presynaptic bouton at any time = rel results in the release of one vesicle. On release, the total mass of glutamate contained within one vesicle exocytosed across the membrane Γpre within the period [rel − ,rel + ], > 0 can be written as

$$\begin{split} M &= \int\_{l\_{\rm rel}-\epsilon}^{l\_{\rm rel}+\epsilon} dt \int\_{\Gamma\_{\rm pre}} D\_e \, \nabla g \cdot \mathbf{n}\_{\Gamma\_{\rm pre}} \, dS \\ &= \int\_{\gamma\_{\rm Vers}-d\_{\rm Vers}/2}^{\gamma\_{\rm Vers}+d\_{\rm Vers}/2} \zeta\_{\rm Vers}^{\rm Gl} \, d\mathbf{y} = \zeta\_{\rm Vers}^{\rm Gl} \, d\_{\rm Vers} \end{split}$$

On average, this total mass is Ves = Ves Gl Ves = 4 3 3 Ves Gl Ves, assuming spherical vesicles with radius Ves = Ves/2 . This yields a mean surface density of

$$
\zeta\_{\rm Vers}^{\rm Gl} = \frac{\pi}{6} c\_{\rm Vers}^{\rm Gl} d\_{\rm Vers}^2 \quad . \tag{5.10}
$$

Released glutamate difuses throughout the synaptic cleft. In accordance with [13], we choose a difusion coefcient of = 3 · 10<sup>5</sup> nm<sup>2</sup> ms−<sup>1</sup> for glutamate in extracellular space. In common reaction kinetic models of synaptic glutamate, the total clearance rate amounts to <sup>c</sup> = 10ms−<sup>1</sup> [15, p.4]. 4/5 of synaptic glutamate is cleared by glial uptake, which leaves only 1/5 being removed by processes on the postsynaptic terminal including uptake and binding to receptors [16]. Furthermore, glutamate gets cleared passively by difusing out of the cleft. We use these proportions to estimate the surface uptake and binding rate densities uptake, bind. Since >> , difusion is the dominant process governing the synaptic glutamate distributions. Thus, we assume → ∞ locally when focusing on the reaction kinetics on the membranes. In this limit, the system is well-mixed and glutamate is homogeneously distributed such that ≡ ∈ R. It behaves like a single-compartment model, which allows us to directly compare the respective fux terms. In particular, we require that all of the clearing fuxes add up to , by

$$\frac{1}{S}k\_{\rm c} \,\mathrm{g}\_{\rm c} = \mathrm{g}\_{\rm c} \int\_{\Gamma\_{\rm post}} k\_{\rm bind} \, dS \tag{5.11}$$

$$\frac{4}{5}k\_{c\text{ }}\mathbf{g}\_{c} = \mathbf{g}\_{c} \int\_{\Gamma\_{\text{a}}} k\_{\text{update}} \, dS \quad . \tag{5.12}$$

We assume that the membrane properties do not change locally, therefore, we assume homogeneous surface rate densities. Thus, we have simple expressions for both surface rate functions,

$$k\_{\rm bind} = \frac{1}{\mathfrak{S}} \frac{k\_c}{h\_l} \tag{5.13}$$

$$k\_{\text{uptake}} = \frac{4}{\mathfrak{F}} \frac{k\_c}{\omega\_s} \quad . \tag{5.14}$$

Here, denotes the width of the synaptic cleft and ℎ describes the full height of the synaptic terminals. Table 5.1 summarizes all parameters and their values associated with our glutamate model.

### **5.2.5 Numerical Methods**

In our validation studies, we solve the glutamate difusion model numerically using the fnite element method with piecewise continuous linear elements for spatial discretization and an implicit Euler scheme for temporal integration.


Table 5.1: Parameters and their values relevant for the glutamate model

Values according to Clements *et al.* [12].

Let <sup>ℎ</sup> denote the set of continuous piecewise linear functions restricted to a discrete triangulation of Ω, such that ∈ <sup>ℎ</sup> vanish at the open boundaries Γopen. Then, the discrete variational form of the glutamate difusion model at time = Δ, where = 0,1, ..., and Δ is the time step, reads: fnd ∈ <sup>ℎ</sup> such that

$$\begin{split} \int\_{\Omega\_{\rm e}} \frac{1}{\Delta t} \mathbf{g} \mathbf{v} + \nabla \mathbf{g} \cdot \nabla \mathbf{v} \, d\mathbf{x} + \int\_{\Gamma\_{\rm a}} k\_{\rm uptake} \left. \mathbf{g} \right|\_{\Gamma\_{\rm a}} \mathbf{v} \, ds + \int\_{\Gamma\_{\rm post}} k\_{\rm bind} \left. \mathbf{g} \right|\_{\Gamma\_{\rm post}} \mathbf{v} \, ds \\ = \int\_{\Omega\_{\rm e}} \frac{1}{\Delta t} \mathbf{g}^{0} \mathbf{v} \, d\mathbf{x} + \sum\_{i} \int\_{\Gamma\_{\rm pre}} v\_{\nu}^{i} \, \delta^{\ast} (t - t\_{i}) \mathbf{v} \, ds \quad \forall \mathbf{v} \in V\_{h} \quad . \end{split} \tag{5.15}$$

Here, 0 refers to the solution at the previous timestep <sup>0</sup> = ( −1)Δ, and ∗ refers to a distribution that, similar to the dirac delta in Equation (5.5), has the property that the glutamate released from each vesicle coincides with its total mass when integrated in time. This distribution may, for example, be a box function which evaluates to 1/ over some interval [, +], or a Gaussian function parameterized by the length scale .

On the other hand, we implement the AMPA model (5.7) using an explicit Euler integration scheme:

$$m\_{\rm AMPA} = m\_{\rm AMPA}^{0} + \Delta t \left( \alpha\_{\rm AMPA} g^{0} \Big|\_{\Gamma\_{\rm psat}} \left( 1 - m\_{\rm AMPA}^{0} \right) - \beta\_{\rm AMPA} m\_{\rm AMPA}^{0} \right) \tag{5.16}$$

The model is implemented and solved using FEniCS [17].

### **5.3 Results**

In the following section, we present the results of our validation study regarding our glutamate model. Figure 5.2 illustrates our experimental setup. It shows the glutamate concentrations at three diferent time points of our numerical solution in the synapse. We start the simulation with the release of a vesicle from Γpre at = 0.01s. In our study, we compare two cases: glutamate release located close to the site of the astrocyte at the bottom of the presynaptic terminal, and release at the opposite site of the terminal, far from the astrocyte. Figure 5.2 displays both cases.

At = 0.01 s, upon release of the glutamate into the synaptic cleft, Figure 5.2 shows highly localized concentrations of glutamate near the regions of release at the presynaptic terminal, indicated by the red regions of high contrast. After 0.1 s, the pronounced glutamate concentration profles have become smoother and the distribution has signifcantly broadened. At = 0.5s, the glutamate distribution has become difusive, and it is clearly visible that glutamate is cleared from the synaptic cleft. It should be noted that the depiction in Figure 5.2 is compressed in height. In fact, the synapse is higher than it is wide by a factor of 20. Note that the apparent increased transport along the horizontal axis is an artifact of this compressed domain. In contrast to our expectations, there seems to be no apparent diference between

Fig. 5.2: Glutamate dynamics in the synaptic cleft. Glutamate distribution in the synaptic cleft is shown at three diferent time points. Vesicles are released at = 0.01 s on the presynaptic terminal (left boundary) either near the astrocyte's location (bottom) or far away from it near the open boundary (top). Glutamate bridges the postsynaptic terminal (to the right) by difusion. For the sake of visualization, the synapse is compressed in height.

the two distributions suggesting that the astrocyte's infuence on the clearance of glutamate might be negligible compared to difusion.

Figure 5.3 summarizes the glutamate fuxes through the membranes of interest and the resulting AMPA-related Na<sup>+</sup> current in a quasi-static approximation of the postsynaptic membrane potential. Figure 5.3A shows the glutamate fux Astro across the astrocyte's membrane Γ<sup>a</sup> as a function of time. After the release near the astrocyte (golden curve), the net fux of glutamate across the membrane increases rapidly and starts to decrease again after 0.25 s as the glutamate gets cleared from the synaptic cleft. In contrast to that, we do not measure a signifcant glutamate uptake through the astrocyte at all if the vesicle is released at the opposite end (teal curve). We note that shortly after the release, Astro becomes negative. According to the boundary condition in Equation (5.3), this corresponds to non-physical negative glutamate concentrations at the boundary. Since ≥ 0, this observation indicates numerical instabilities that may be attributed to the Gibbs phenomenon. Despite that, the concentration and Astro quickly stabilize and exhibit the expected behaviour.

Fig. 5.3: Flux dynamics vs. time in the glutamate model. **a)**: Total glutamate uptake by the astrocyte. **b)**: Total glutamate fux through the postsynaptic terminal, Post, due to binding to the AMPA receptor. **c)**: Total Na<sup>+</sup> current AMPA as a result of AMPA activation caused by the arriving glutamate in case of a quasi-static membrane potential at post = −65mV.

Figure 5.3B depicts the glutamate fux through the postsynaptic membrane, Post, in time. Post has important implications for the excitation of the postsynaptic neuron because it indirectly refects the activation of the AMPA receptors that control Na<sup>+</sup> channels on the membrane. The profle of glutamate uptake through the postsynaptic terminal is qualitatively similar to that of the astrocyte. On the falling branch of the curve however, we see that if the astrocyte is locally present, less glutamate reaches the postsynaptic terminal as a result of active clearance by the astrocyte. This validates our model with respect to the literature [3, 5, 4].

Figure 5.3C shows the magnitude of the AMPA-regulated Na<sup>+</sup> current across the postsynaptic terminal in time. Thus, it refects the total activity of the AMPA receptors on Γpost as glutamate approaches that membrane and binds to the receptors. AMPA is calculated from the AMPA gating variable, AMPA, using Equation (5.8). We assume that in the simulated period, the postsynaptic membrane potential has not deviated too far from its resting potential despite the fnite AMPA current. With this quasi-stationary approximation, we justify Post ≃ −65mV, the neuron resting potential, and are able to give a coarse estimate of AMPA without having to couple our model. Similar to Tewari *et al.* [14], we choose AMPA = 0.35nS and AMPA = 0mV. Negative signs indicate a net Na<sup>+</sup> current out of the dendrite instantiating depolarization across the membrane. As glutamate reaches the postsynaptic neuron, we see a growing increase in the AMPA-specifc Na<sup>+</sup> current, as expected. In contrast to Post, AMPA does not saturate within the simulated time window. We note that while less glutamate reaches the postsynaptic neuron in the region of infuence of the astrocyte (Figure 5.3B), unexpectedly, the resulting AMPA current AMPA of this region appears to be more pronounced.

### **5.4 Discussion**

We developed a difusion-based, spatio-temporal model of the neurotransmitter glutamate in glutamatergic hippocampal synapses. We carried out validation studies demonstrating that our model reproduces the astrocyte's role in glutamate uptake in the tripartite synapse. In addition, we observed that astrocytic glutamate uptake translates into a decrease in glutamate fux across the postsynaptic terminal, showing that astrocytes temper the excitatory signal that is transmitted to postsynaptic neurons. This is a widely recognized key role of astrocytes [3, 5, 4]. We also saw in our model that the binding of glutamate to AMPA receptors in the postsynaptic neuron results in a net AMPA-specifc Na<sup>+</sup> current AMPA, initiating depolarization in the postsynaptic neuron. Thus, our model successfully translates synaptic chemical cues into an electrical signal in the postsynaptic neuron. Unexpectedly, we observed that the total postsynaptic glutamate fux, Post, inversely afects the strength of the AMPA-specifc Na<sup>+</sup> current. This efect could be explained by the astrocyte causing a narrower glutamate distribution, resulting in a locally restricted but higher AMPA activation. As a consequence, the total Na<sup>+</sup> current might be stronger. The glutamate profles in Figure 5.2, however, indicate otherwise. Moreover, our quasi-static approximation of AMPA does not resemble the actual Na<sup>+</sup> current dynamics accurately. Using the correct dynamic potentials, e.g. by coupling our glutamate model to the KNP-EMI model, might yield diferent results. Finally, we observed characteristics in our simulation results that indicate issues of numerical stability near the astrocyte during release. Similar instabilities might be causing the discrepancy in AMPA that we can not explain in the scope of our model. Thus, the next logical step would be to study the stability and convergence of our methods comprehensively.

Furthermore, we observed a signifcant diference in postsynaptic activity with and without the astrocyte. However, the efect is not as pronounced as we would have expected. Conversely, it appears rather minute. The reason for that might lie in our choice of the astrocyte's location in terms of distance from the neurons. In our model, we have chosen a distance of 50 nm, which is 2.5-fold larger than the width of the synaptic cleft. As a next step, we suggest situating the astrocyte closer to the neurons. Lastly, the astrocyte's role of clearing the synaptic cleft is expected to be most signifcant when it is saturated with glutamate, for example, in bursting neurons. In this case, excessive concentrations of neurotransmitters can have a neurotoxic efect. In our study, we only considered single-vesicle release, but our model accounts for an arbitrary amount. Consequently, with longer simulated time windows, multi-vesicular release studies resembling spike trains of arbitrary frequency are an obvious and promising connecting factor for future work.

Other computational studies measuring glutamate fux in a synapse could not be found in the literature. However, some physiological studies may be used to contextualize our results. Although glutamate fux across an astrocyte is not directly measured *in vitro*, it can be carried out indirectly by measuring the changing glutamate concentration of a medium containing astrocyte cultures. In this regard, the studies by Farinelli and Nicklas [18] and Fonseca *et al.* [19] show similar glutamate clearance dynamics as our results (e.g., Figures 5.2 and 5.3, respectively). These astrocytic clearance dynamics are measured over a much longer time period. However, the observed pattern may be extended to our study in that there is an initial fast rate of glutamate clearance, followed by a plateau stage during which the rate of clearance slows as the relative glutamate concentration approaches 0 mM.

A clear extension to our study would be utilizing the model in tandem with the KNP-EMI model. Our model design integrates well into EMI or KNP-EMI frameworks by adding the AMPA-specifc contribution to the Na<sup>+</sup> current, AMPA. In the EMI and KNP-EMI frameworks (e.g. [8]), the choice of the functional form of ion is subject to modelling, usually in form of a solution to a Hodgkin-Huxley type ODE system [20]. When integrating our model into the KNP-EMI framework, we suggest setting ion = AMPA on the postsynaptic terminal (Γpost in Figure 5.1) and Hodgkin-Huxley based models on the remaining postsynaptic membranes (Γd).

Considering further astrocyte-related mechanisms, such as the dual function of glutamate uptake and release [6], the KNP-EMI model is expected to be especially well-suited for modelling glutamate uptake. As a frst step into modeling neurotransmitters within the KNP-EMI framework, we designed our model with simplicity in mind. For this reason, our model only features what we identifed as the most relevant processes for signal transmission, excluding presynaptic glutamate re-uptake and astrocytic glutamate release. These mechanisms should, however, be considered in a refned model, because they are supposedly relevant for realizing glutamate bufers in the synapse. Additionally, the dual function mechanism is controlled in part by intracellular pH levels [21], which may be well handled by the ionic concentration and difusion components of the KNP-EMI model.

In summary, we successfully developed a glutamate extension to the KNP-EMI model in the tripartite synapse. Our model exhibits all the relevant mechanisms necessary for signal transmission between neurons across the synaptic cleft. In our model, the astrocyte modulates the signaling molecule concentrations as expected. Their efect could be more pronounced, but we expect that our choice regarding the astrocyte's distance from the synapse had a signifcant infuence on this outcome. Additionally, the resulting Na<sup>+</sup> current behaves contrarily to expectations in response to the astrocyte's infuence. Future work on our model should address our concerns regarding this particular parameter choice and the numerical stability of the model.

5 EMI for Astrocyte-Neuron Interactions 77

### **References**


21. Julia Langer, Niklas J Gerkau, Amin Derouiche, Christian Kleinhans, Behrouz Moshref-Ravasdjani, Michaela Fredrich, Karl W Kaftz, Gerald Seifert, Christian Steinhauser, and ¨ Christine R Rose. Rapid sodium signaling couples glutamate uptake to breakdown of atp in perivascular astrocyte endfeet. *Glia*, 65(2):293–308, 2017.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 6 Inducing Flow Instabilities in Aneurysm Geometries via the Reynolds-Orr Method**

Alessandro Contri, Christina Taylor, Justin Tso, and Ingeborg Gjerde

**Abstract** Ruptured intracranial aneurysms are the leading cause of hemorrhagic strokes. Although intracranial aneurysms are prevalent in as much as 4% of the adult population, most aneurysms do not rupture, and their growth and risk of rupture remains difcult to predict. Previous studies have identifed abnormal wall shear stress patterns and blood pressure spikes, both features of unstable fows, as key factors in aneurysm behavior. While computational fuid dynamics has been enlisted to help study risk of rupture, no consensus has been established on how to impose unstable fow conditions. Here, we use Reynolds-Orr instability analysis to calculate fow perturbations that are capable of inducing unstable fow in 3D arterial geometries with and without aneurysms. These perturbations lay the foundation for future wall shear stress studies by providing mathematically consistent conditions for time-dependent Navier-Stokes simulations.

### **6.1 Introduction**

Intracranial aneurysms are relatively common in adults, with a global prevalence of 2-3%. Intracranial aneurysm rupture leads to aneurysmal subarachnoid hemorrhage, which is fatal in about 60% of cases and carries major risk of brain damage and

Alessandro Contri

Department of Mathematical Sciences, Norwegian University of Science and Technology, NO

Christina Taylor

Justin Tso

Department of Biomedical Engineering, Johns Hopkins University, USA

Ingeborg Gjerde (corresponding author)

Department of Computational and Applied Mathematics and Operations Research, Rice University, USA

Department of Numerical Analysis and Scientifc Computing, Simula Research Laboratory, NO e-mail: ingeborg@simula.no

disability for survivors [1]. However, only 1-2% of intracranial aneurysms rupture annually, and most are asymptomatic [2]. The decision to treat unruptured aneurysms must therefore be carefully weighed against the risk of rupture, a metric that remains difcult to quantify.

Aneurysm growth and rupture have been linked to local hemodynamic conditions, with wall shear stress (WSS) primary among them [3]. Computational fuid dynamics (CFD) studies allow for the high-resolution investigation of WSS in patient-specifc arterial geometries through simulations of the Navier-Stokes equations. Several CFD studies have investigated the role of WSS in aneurysm rupture with conficting results. While some studies have attributed risk of rupture to high WSS at thin and hard aneurysm walls, more recent studies have found strong relationships with low WSS at aneurysm rupture points (reviewed in [4]). Meng *et al.* have noted that both high and low WSS may increase rupture risk: high WSS causes degeneration of the cellular matrix and cell apoptosis, while low WSS and recirculation promotes an infammatory environment that weakens the arterial wall and drives aneurysm growth [5].

Several studies have found rapid fuctuations in blood pressure to be a common risk factor for aneurysm rupture. Wouter *et al.* found that 43% of ruptures occurred during stressful events, while Vlak *et al.* found high population-attributable risks for cofee drinking and vigorous physical exercise [6, 7]. Matsuda *et al.* also found high incidence of rupture in activities associated with the Valsalva maneuver, which results in sudden pressure changes across the aneurysm wall [8]. These perturbations in fow are highly relevant to aneurysm rupture, but are not accounted for in most CFD studies.

These fow instabilities can be investigated using defnitions of instability based on kinetic energy, which originated with Reynolds and Orr [9]. Analysis using these defnitions, as revisited by Scott, provides an exact relation between a base fow and the evolution of the kinetic energy of a perturbation. The method is nonlinear with no approximations and yields a well-posed linear symmetric eigenvalue problem that can be solved with standard methods. In the resulting eigenpairs the fow perturbations are captured in the eigenvectors, while negative eigenvalues indicate the potential for instability with respect to kinetic energy. In unstable eigenpairs, the magnitude of the eigenvalues indicates the initial growth rates of the perturbed kinetic energy. However, the most energetically unstable eigenpair is a function of both the magnitude of its eigenvalue and its eigenvector [9].

While the Reynolds-Orr eigenpairs may not exactly refect physiological fow instabilities, they provide a reproducible way to induce fow instability. The readily available mathematical framework behind this approach allows this analysis to easily be reproduced in diferent arterial geometries. In this work, we investigate some of the nuances of using this approach to induce unstable fow on arterial geometries with and without aneurysms taken from the Aneurisk data set from Emory University [10].

In the remainder of this paper, we will frst discuss the Reynolds-Orr method in more detail in Section 6.2.1 followed by the details of our time-dependent Navier-Stokes solver in Section 6.2.2. Next, we will present some of our initial fndings on using the Reynolds-Orr to induce instability in Section 6.3 before concluding with a discussion on our fndings in Section 6.4.

### **6.2 Methods**

In this section we consider the mathematical equations and the subsequent implementation of the problem at hand, part of which was previously presented in [11]. The artery is considered as a pipe with a no-slip wall and prescribed pressures at the inlet and outlets. The absolute pressure is disregarded as it does not change the relative pressure drop in the model.

### **6.2.1 Reynolds-Orr Instability**

The derivation of the eigenvalue problem from which we hope to detect the most unstable modes, given the base fow, follows the one presented in [9]. The derivation is nontrivial for pressure boundary conditions; therefore, we consider here only Dirichlet boundary conditions for the velocity. Let (**u**, ) be the solution to the following unsteady Navier-Stokes equations:

$$
\partial\_t \mathbf{u} - \nu \Delta \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p = \mathbf{f} \quad \text{in} \quad \Omega \times (0, T], \tag{6.1}
$$

$$
\nabla \cdot \mathbf{u} = 0 \quad \text{in} \quad \Omega \times (0, T], \tag{6.2}
$$

$$\mathbf{u} = \mathbf{g} \quad \text{on} \quad \partial \Omega \times (0, T], \tag{6.3}$$

$$\mathbf{u} = \mathbf{u}^0 \quad \text{on} \quad \Omega \times \{0\}. \tag{6.4}$$

Let now the couple (**w**, ) solve the same equations (6.4) with the same boundary condition **g**, but diferent initial data **w** 0 , such that **w** <sup>0</sup> ≠ **u** 0 . Defne (**v**, ) as the diference between the two solutions:

$$\mathbf{v} = \mathbf{u} - \mathbf{w}, \qquad o = p - q \,. \tag{6.5}$$

Then, taking the two versions of (6.4) with the diferent initial data and subtracting them, (**v**, ) satisfy:

$$
\partial\_t \mathbf{v} - \nu \Delta \mathbf{v} + (\mathbf{u} \cdot \nabla \mathbf{u} - \mathbf{w} \cdot \nabla \mathbf{w}) + \nabla o = \mathbf{0} \quad \text{in} \quad \Omega \times (0, T], \tag{6.6}
$$

$$
\nabla \cdot \mathbf{v} = 0 \quad \text{in} \quad \Omega \times (0, T], \qquad (6.7)
$$

$$\mathbf{v} = \mathbf{0} \quad \text{on} \quad \partial \Omega \times (0, T], \qquad (6.8)$$

$$\mathbf{v} = \mathbf{v}^0 = \mathbf{u}^0 - \mathbf{w}^0 \quad \text{on} \quad \Omega \times \{0\}. \qquad (6.9)$$

Following [9], the nonlinear term is rewritten as

$$\mathbf{u} \cdot \nabla \mathbf{u} - \mathbf{w} \cdot \nabla \mathbf{w} = \mathbf{u} \cdot \nabla \mathbf{u} - \mathbf{u} \cdot \nabla \mathbf{w} + \mathbf{u} \cdot \nabla \mathbf{w} \mathbf{w} \cdot \nabla \mathbf{w} \tag{6.10}$$

$$\mathbf{v} = \mathbf{u} \cdot \nabla \mathbf{v} + \mathbf{v} \cdot \nabla \mathbf{w} = \mathbf{u} \cdot \nabla \mathbf{v} + \mathbf{v} \cdot \nabla \mathbf{u} - \mathbf{v} \cdot \nabla \mathbf{v} \tag{6.11}$$

Multiplying (6.9) by **v**, using (6.11) and integrating over Ω yields

$$a\left(\partial\_t \mathbf{v}, \mathbf{v}\right) + a\left(\mathbf{v}, \mathbf{v}\right) + c\left(\mathbf{u}, \mathbf{v}, \mathbf{v}\right) + c\left(\mathbf{v}, \mathbf{u}, \mathbf{v}\right) - c\left(\mathbf{v}, \mathbf{v}, \mathbf{v}\right) - b\left(\mathbf{v}, o\right) = 0. \tag{6.12}$$

where (·, ·) denotes the standard 2 (Ω) inner product and

$$a(\mathbf{u}, \mathbf{v}) = \nu \int\_{\Omega} \nabla \mathbf{u} : \nabla \mathbf{v}, \tag{6.13}$$

$$c(\mathbf{u}, \mathbf{v}, \mathbf{w}) = \int\_{\Omega} (\mathbf{u} \cdot \nabla \mathbf{v}) \cdot \mathbf{w},\tag{6.14}$$

$$b(\mathbf{u}, p) = \int\_{\Omega} \nabla \cdot \mathbf{u} p. \tag{6.15}$$

Defning now the following space:

$$\mathbf{V} \coloneqq \{ \mathbf{v} \in H^1(\mathfrak{Q})^3 : \nabla \cdot \mathbf{v} = 0 \quad \text{and} \quad \mathbf{v} = 0 \quad \text{on} \quad \partial \Omega \} \tag{6.16}$$

it can be observed from [9] that for any **w** ∈ **V** one have

$$c\left(\mathbf{w}, \mathbf{v}, \mathbf{v}\right) = 0, \quad b\left(\mathbf{v}, o\right) = 0. \tag{6.17}$$

(6.12) can then be rewritten in the form

$$\frac{1}{2}\frac{\partial}{\partial t}\int\_{\Omega}|\mathbf{v}|^{2} = -\nu \int\_{\Omega} |\nabla \mathbf{v}|^{2} - \int\_{\Omega} (\mathbf{v} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} \tag{6.18}$$

$$\mathbf{u} = -\nu \int\_{\Omega} |\nabla \mathbf{v}|^2 - \frac{1}{2} \int\_{\Omega} \mathbf{v}^T (\nabla \mathbf{u} + \nabla \mathbf{u}^T) \mathbf{v}. \tag{6.19}$$

The frst term of the equation can be seen as the time derivative of the kinetic energy of the fow and the right hand side describes its evolution. Following standard arguments (i.e. integration in time of the left hand side etc.) we can say that the fow **u** is energy unstable at time = 0 if there exists a **v**<sup>0</sup> ∈ **V** such that

$$-\frac{1}{2}\int\_{\Omega} \mathbf{v}\_0^T (\nabla \mathbf{u}\_0 + \nabla \mathbf{u}\_0^T) \mathbf{v}\_0 - \nu \int\_{\Omega} |\nabla \mathbf{v}\_0|^2 > 0. \tag{6.20}$$

At this point we defne

$$\lambda\_{\mathbf{v}} = \frac{B\_{\mathbf{u}}(\mathbf{v}, \mathbf{v})}{a(\mathbf{v}, \mathbf{v})} \quad \text{with} \quad B\_{\mathbf{u}}(\mathbf{v}, \mathbf{w}) = \frac{1}{2} \int\_{\Omega} \mathbf{v}^{T} (\nabla \mathbf{u} + \nabla \mathbf{u}^{T}) \mathbf{w} \tag{6.21}$$

Comparing (6.19) with (6.21) we can reformulate an equivalent instability condition. Namely, the fow is unstable at = 0 if there is a **v**<sup>0</sup> ∈ **V** such that

$$
\lambda\_{\mathbf{v}\_0} < 1.\tag{6.22}
$$

Given that (**v**,**v**) = ∥∇**v**∥ 2 , we have in particular that if ∃**v** such that = inf0≠**v**∈**<sup>V</sup> <sup>v</sup>** < −1 and thus

$$
\int \left(\frac{1}{\|\nabla \mathbf{v}\|^2}\right) \frac{\partial}{\partial t} \int\_{\Omega} |\mathbf{v}|^2 = -2\nu(1 + \lambda\_{\mathbf{v}}) > 0,\tag{6.23}
$$

the fow is unstable. Recalling Poincare's inequality, stating ´ ∥**v**∥ <sup>2</sup> ≤ ∥∇**v**∥ <sup>2</sup> ∀**v** ∈ **V**, we can deduce

$$\frac{\partial}{\partial t} \log \left( \int\_{\Omega} |\mathbf{v}|^2 \right) \ge -\frac{2\nu}{C\_P} (1 + \lambda\_\mathbf{v}).\tag{6.24}$$

We conclude that the most negative **<sup>v</sup>** leads to the most unstable mode. However, as noted in [9], the most unstable mode may not be the most persistent. Additionally, [9] proves that the solution to = inf0≠**v**∈**<sup>V</sup> <sup>v</sup>** solves the eigenvalue problem (,**v**) ∈ (R,**V**) such that

$$B\_{\mathbf{u}}(\mathbf{v}, \mathbf{w}) = \lambda a(\mathbf{v}, \mathbf{w}) \quad \forall \mathbf{w} \in \mathbf{V}. \tag{6.25}$$

The detailed derivation can be found in [9].

### **6.2.2 Unsteady Navier-Stokes FEM Discretization**

We looked to study the evolution in time of the potentially unsteady modes detected using the method described previously. In order to progress in time we needed to implement a time-dependent Navier-Stokes solver. To do so, we used a splitting method, which despite its low accuracy in time is fast and easy to implement from the computational point of view. The motivation behind splitting schemes is to consider the frst two equations in 6.4 separately. We chose a modifed version of Chorin's method [12], the so-called incremental pressure correction scheme (IPCS) [13], which gives improved accuracy compared to the original scheme at little extra cost.

The IPCS scheme involves three steps. First, we compute a tentative velocity **u** ∗ by advancing the momentum equation (6.4) by a midpoint fnite diference scheme in time, using the pressure from the previous time interval. We also linearize the nonlinear convective term by using the known velocity **u** from the previous time step: **u** · ∇**u** . The variational problem for this frst step is

$$((\mathbf{u}^\* - \mathbf{u}^n)\Delta t, \mathbf{v}) + (\mathbf{u}^n \cdot \nabla \mathbf{u}^n, \mathbf{v}) + a(\mathbf{u}^{n + \frac{1}{2}}, \mathbf{v}) - b(\mathbf{v}, p^n) \tag{6.26}$$

$$+(p^n\mathbf{n},\mathbf{v})\_{\partial\Omega} - \nu(\nabla \mathbf{u}^{n+\frac{1}{2}} \cdot \mathbf{n}, \mathbf{v})\_{\partial\Omega} = (\mathbf{f}^{n+1}, \mathbf{v})\tag{6.27}$$

where

$$\mathbf{u}(\mathbf{u}, \mathbf{v})\_{\partial\Omega} = \int\_{\partial\Omega} \mathbf{u} \mathbf{v} \, \mathrm{d}s, \quad \mathbf{u}^{n+\frac{1}{2}} = (\mathbf{u}^n + \mathbf{u}^{n+1})/2. \tag{6.28}$$

Since we are modeling fow in a "pipe" (the artery) with known infow and outfow, the scheme has the hidden assumption that the derivative of the velocity in the direction of the channel is zero at the infow and outfow, corresponding to a fow that is "fully developed," or that the fow feld doesn't change signifcantly upstream or downstream of the domain.

The second step is to use the computed tentative velocity to compute the new pressure :

$$(\nabla p^{n+1}, \nabla q) = (\nabla p^n, \nabla q) - \Delta t^{-1} (\nabla \cdot \mathbf{u}^\*, q) \tag{6.29}$$

:

Finally, we compute the corrected velocity **u** +1

$$\mathbf{u}(\mathbf{u}^{n+1}, \mathbf{v}) = (\mathbf{u}^\*, \mathbf{v}) - \Delta t(\nabla(p^{n+1} - p^n), \mathbf{v}).\tag{6.30}$$

In summary, we solve the incompressible Navier–Stokes equations efciently by solving a sequence of three linear variational problems in each time step. The FEM space chosen to solve the problem is (**u**ℎ, ℎ) ∈ (**V** 2 <sup>Γ</sup>,ℎ,**<sup>V</sup>** 1 ℎ ) where, if we let **V** ℎ denote the space of piecewise polynomials of degree on a regular mesh of the domain Ω, the vector valued polynomial space is

$$\mathbf{V}\_{\Gamma,h}^k = \{ \mathbf{v} \in (\mathbf{V}\_h^k)^3 : \mathbf{v}|\_{\partial \Omega \mid \Gamma} = 0 \}. \tag{6.31}$$

### **6.3 Results**

For our numerical experiments we used a single artery mesh both with and without a aneurysm. These geometries are shown in Figure 6.1.

Fig. 6.1: The artery geometry without (Case 0) and with (Case 1) an aneurysm used for our numerical experiments.

As previously stated, eigenvalues with negative sign and magnitude greater than one are necessary to induce energy instability. Eigenpairs were computed using a shifted power iteration with LU preconditioning to target eigenvalues near -1. We observed that without a sufcient pressure drop, eigenvalues less than or equal to -1 were not guaranteed. Using Δ = 600Pa (4.5mmHg) and 1000Pa (7.5mmHg) both yielded at least two eigenvalues less than -1 in both the healthy artery and terminal aneurysm. Tables 6.1 and 6.2 show the eigenvalues reported by the solver. Figures 6.2 and 6.3 show example perturbations associated with the Δ = 600Pa pressure drop case.


Table 6.1: Eigenvalues for the healthy artery without an aneurysm for diferent pressure drop values.

Table 6.2: Eigenvalues for the artery with an aneurysm for diferent pressure drop values.


<sup>∗</sup> The 600Pa pressure drop did not yield a sixth eigenvalue near ±1.

With the eigenpairs in hand, it remained to simulate the kinetic energy evolution of the eigenpairs with eigenvalues less than -1. The Navier-Stokes equations are notoriously difcult to simulate especially in the face of physical instabilities. In the case of our scheme a severely prohibitive time step was needed to ensure numerical stability. The prohibitive time step coupled with a lack of computing resources

Fig. 6.2: Example of the base fow without and with perturbation added in the healthy artery case for Δ = 600Pa.

Fig. 6.3: Example of base fow without and with perturbation added in the aneurysm case for Δ = 600Pa.

restricted us from being able to perform in-depth analysis of the kinetic energy evolution. However, we were able to calculate the kinetic energy evolution for the above eigenpairs over a very short time period. Figures 6.4 and 6.5 show the L2 norm of the velocity (the kinetic energy) of the solutions over a time period of 0.1 seconds.

### **6.4 Discussion**

From our initial results, only four eigenpairs seem to induce unstable fow as indicated by increasing kinetic energy despite there being 16 negative eigenvalues. However, as shown in Figure 6.4 for the 1000Pa case with = −1.319, it seems possible for the kinetic energy to initially decrease before increasing. This is not completely unexpected as [9] showed the kinetic energy can oscillate, though in their experiments

Fig. 6.4: Evolution of the L2 norm of the velocity in time for healthy artery eigenpairs.

this oscillation was coupled to an overall decay after an initial increase in energy. Our short simulation time may not have been sufcient to capture instabilities that require more time to develop. Additionally, as was noted in [9] a larger magnitude negative eigenvalue does not necessarily guarantee the most unstable mode due to the instability's additional dependence on the gradient of the perturbed velocity. In Figure 6.5 this can be seen in the 1000Pa case by the eigenpair with = −0.785 inducing an initial increase in energy while the eigenpair with = −1.154 did not.

While the restrictive time step required by the Navier-Stokes solver did not allow us to complete all initial goals of this project, we have been able to show some preliminary results showing that perturbations computed using the Reynolds-Orr method can induce energy instability for a given geometry. Future work would entail longer simulations to confrm whether all negative eigenvalues with sufcient

Fig. 6.5: Evolution of the L2 norm of the velocity in time for the aneurysm eigenpairs.

magnitude do lead to instability eventually. Additionally, we have completed the work started in [11] by correcting errors that had previously made when calculated the eigenmodes. The codes produced here provide the framework needed to continue studying aneurysm rupture risk through the Reynolds-Orr instability analysis when provided sufcient computing resources.

6 Inducing Flow Instabilities in Aneurysms 89

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 7 Impact of Pathological Vascular Remodelling on Right Ventricular Mechanics**

Jiaxin (Katie) Cui\*, Mariluz Rojo Domingo\*, Ryan Konno\*, Claudia A Manetti\*, George Kagugube\*, Oscar Odeigah, and Joakim Sundnes \* These authors contributed equally to this work

**Abstract** Pulmonary arterial hypertension (PAH) is a rare disorder characterized by elevated blood pressure and pulmonary vascular resistance, often followed by right ventricular hypertrophy and heart failure. The efect of PAH and its treatments on the mechanics, function, and remodelling of the right ventricle (RV) is currently not well understood. To study cardiac biomechanics and functionality as PAH progresses, we implemented a computational model of the heart simulating right ventricular maladaptive remodelling. Our Windkessel-based model, which accounts for direct ventricular interaction and the presence of the pericardium, is utilized to simulate various disease stages of PAH. We fnd that the pericardium has a larger efect on heart performance than ventricular interaction through the septum. We also examined the efectiveness of two treatments, ventricular assist device (RVAD) and atrial septostomy, on diseased hearts. We show that while both pulsatile and continuous RVADs restore cardiac function, pulsatile RVAD improves cardiac output 29.4% more than continuous RVAD. We also demonstrate that atrial septostomy improves

Mariluz Rojo Domingo Department of Bioengineering, University of California, San Diego, USA

Ryan Konno Department of Mathematics, Simon Fraser University, CA

Claudia A Manetti Maastricht University, NL

George Kagugube University College London, UK

Oscar Odeigah Department of Computational Physiology, Simula Research Laboratory, NO

Joakim Sundnes (corresponding author) Department of Computational Physiology, Simula Research Laboratory, NO e-mail: sundnes@ simula.no

Jiaxin (Katie) Cui

Department of Biomedical Engineering, Viterbi School of Engineering, University of Southern California, USA

cardiac output by 19.5%. Our model can be further extended by simulating the heart's response to other treatments such as extracorporeal membrane oxygenation (ECMO), and by incorporating ventricular remodelling growth simulations and fnite-element ventricular modelling.

### **7.1 Introduction**

Pulmonary arterial hypertension (PAH) is characterized by high blood pressure and high pulmonary vascular resistance in the arteries in the lungs [1]. As the pressure and resistance in the pulmonary arteries and capillaries rises, the right ventricular afterload is elevated. Over time, there is a larger demand for work done by the heart to achieve the same level of cardiac output as that of a healthy one. To increase the amount of contractile force generated by the heart, the right ventricle (RV) undergoes abnormal enlargement and increases its wall thickness, stifness, and contractility. This phenomenon, which is observed in almost all PAH patients, is referred to as right ventricular hypertrophy [2]. Moreover, the heart compensates for the elevation in pulmonary arterial pressure and resistance via vascular proliferation and remodelling of small pulmonary arteries [3]. In advanced stages of PAH, the severe dilation of the RV causes septal bowing, i.e., a leftward shift of the septum that lowers left ventricular volume and limits left ventricular flling. Hence, left ventricular output is reduced and, in turn, cardiac output and mean arterial pressure decrease [4]. All these compensatory mechanisms enable the heart to maintain its performance in the early stages of PAH; however, sustained structural and functional changes in the RV can eventually become detrimental to cardiac function and may even result in right heart failure and death of a patient [5].

Although life-threatening, PAH remains incurable. Treatment strategies are available to alleviate the vascular symptoms (e.g., vasoconstriction) [6], and to delay disease progression by unloading the RV. Therapies proven to be efective in PAH patients include: drugs in the form of vasodilators, inotropes and anticoagulants; extracorporeal membrane oxygenation (ECMO) devices; mechanical support via a right ventricular assist device (RVAD); atrial septostomy; or lung transplantation. [7]. Multiple *in situ* and *in silico* investigations focus on PAH, yet there is still a gap in our understanding of the biomechanic and hemodynamic efect of PAH on the RV and its subsequent remodeling (and failure). [3, 4, 8, 9, 10] Furthermore, the progressive remodeling of the RV is inevitably infuenced by its surrounding environment; in other words, the septum, the pericardium, the left ventricle (LV), and barorefex must be considered in order to investigate RV performance in PAH as accurately as possible. Some computational studies have already looked into changes in heart function with ventricular interdependence under the infuence of pericardium and barorefex. [11, 12, 13] Nonetheless, most rigorous biomechanical and hemodynamic models have yet to be utilized for research involving treatments of PAH and right ventricular remodelling. Our objective was to develop a realistic computational model that simulates the heart conditions in a PAH setting to gain insight into the various maladaptations characteristic of this condition as well as the efect of PAH treatments on RV mechanics. Our theoretical analysis supplements experimental research techniques to further enhance our knowledge of PAH therapies and how the progressive adverse remodeling impacts the biomechanics and hemodynamics of the heart.

### **7.2 Methods**

In this work, we implement a closed-loop, lumped-parameter model of the heart and the circulatory system [8, 12] to study the impact of increased pulmonary vascular resistance on RV hemodynamics. This Windkessel-based model captures the efects of PAH at diferent severities while accounting for direct ventricular interaction via the interventricular septum, the pericardium constraint, barorefex, and the heart's dynamic growth and remodelling in response to varying loads. Moreover, our model is utilized to examine two treatment options for PAH-induced RV dysfunction. In particular, we investigate the impact of RVAD and atrial septostomy on cardiovascular hemodynamics in the PAH setting.

### **7.2.1 Base Model of PAH**

Our model of the cardiovascular system combines elements of the cardiac anatomy with those of the systemic and pulmonary circulations. The pressure, volume, and blood fow through the circulatory system are represented using a set of ordinary diferential equations. The relationship between the pressure and volume for each component in the circulatory system can be represented linearly by

$$V = CP + V\_{ds}$$

where is volume, is pressure, is the compliance of the vessel, and is the volume of the vessel at zero pressure. While and of the arteries and veins are typically constant (unless there are physiological or pathological changes), and of the heart chambers are time-varying periodic functions that represent active contraction [14]. The circulatory system model consists of individual elastic chambers and a series of resistors and diodes with fxed and variable capacitors branching from the main circuit (Figure 7.1) [15]. The systemic and pulmonary vasculatures are divided into their corresponding resistive and fxed capacitive components due to their elastic nature. Heart valves reinforce uni-directional blood fow, so they are represented as diodes. In addition, the mitral and tricuspid valves exhibit non-negligible resistive elements to block blood fow between heart chambers. Systole and diastole, the two phases of the cardiac cycle, are captured by variable capacitors that represent the heart chambers.

Fig. 7.1: Base Windkessel circulation model describing the coupling between the heart and the systemic and pulmonary circulation, including RVAD and atrial septostomy treatments.

The volume of each individual compartment is related to the rate of infow and outfow via

$$\frac{dV\_i}{dt} = \mathcal{Q}\_{in} - \mathcal{Q}\_{out},\tag{7.1}$$

where is the volume of chamber , and and represent the fow in and out of the chamber, respectively. We approximate the blood fow at each component, which is driven by pressure diference, to be linearly dependent on pressure through the following relation:

$$\mathcal{Q} = \frac{P\_{upstream} - P\_{downstream}}{R},\tag{7.2}$$

where is the resistance of the given connection.

Ventricular and atrial pumping characteristics are modelled by modifying timevarying elastance, (). Mathematically, instantaneous ventricular pressure () can be related to instantaneous ventricular volume () via

$$P(t) = \left[P\_{es}(V) - P\_{ed}(V)\right]e(t) + P\_{ed}(V),\tag{7.3}$$

where

$$P\_{es}(V) = E\_{es}(V - V\_0),\tag{7.4}$$

$$P\_{ed}(V) = \alpha \left( e^{\beta(V - V\_0)} - 1 \right),\tag{7.5}$$

7 Right Ventricular Mechanics 95

$$e(t) = \begin{cases} \frac{1}{2} \left[ \sin \left( \frac{\pi}{T\_{ex}} t - \frac{\pi}{2} \right) + 1 \right], & \text{if } 0 < t < \frac{3T\_{ex}}{2} \\\\ 0.5 \exp \left\{ \left( t - \frac{3T\_{ex}}{2} \right) / \tau \right\}, & \text{otherwise.} \end{cases} \tag{7.6}$$

Here, is end-systolic pressure, is end-diastolic pressure, is end-systolic elastance, is time required to achieve end-systole, is the time constant of relaxation, and 0, , , are constant model parameters in the end-diastolic pressurevolume relationship. The parameters of the base model were adapted, based on Punnoose *et al.* [8], to capture diferent stages of PAH (mild, moderate, severe, and very severe, i.e., cardiogenic shock). We initialized the model with standard parameters for a healthy heart, and then varied its parameters to match those of the corresponding phase. When the parameters were altered for each diseased state, common RV function metrics, such as the stroke volume, the ejection fraction and cardiac output, changed too. We could thus observe a distinct pressure-volume relationship for each degree of right-sided heart failure before implementing any treatments. Overall, the most relevant parameters in the circulation model include systemic and pulmonary vessel resistance, compliance and elastance. The varying efects of most of these parameters on the model's behaviour was further explored by performing a sensitivity analysis.

### **7.2.2 Ventricular Interaction and the Pericardium**

In this study, we aimed to improve the realism of the previously described lumpedparameter model by accounting for the interventricular interaction via the septum wall. The left and the right ventricles act independently of each other; their timevarying volume-elastance relation is represented by the system shown above in Equation 7.6. Following Santamore *et al.* [11], three time-varying elastances (right ventricular free wall , left ventricular free wall , and septum wall ) were considered, so that the interdependence between the two ventricles was simulated by our model. The pressures of both ventricles at the end of diastole () becomes

$$P\_{edLV}(V) = \alpha \left( e^{B\_{lv}(V - V\_0)} - 1 \right),\tag{7.7}$$

where

$$B\_{l\boldsymbol{\nu}} = B\_{l0} + m\_{r\boldsymbol{l}} V\_{r\boldsymbol{\nu}},\tag{7.8}$$

and

$$P\_{edRV}(V) = \alpha (e^{B\_{rv}(V - V\_0)} - 1),\tag{7.9}$$

where

96 Right Ventricular Mechanics

$$B\_{r\upsilon} = B\_{r0} + m\_{lr} V\_{l\upsilon} \tag{7.10}$$

<sup>0</sup> and <sup>0</sup> represent the ventricular coefcients when the volume is equal to zero. is the sensitivity of to changes in right ventricular volume while is the sensitivity of to changes in left ventricular volume. The pressure development in each ventricle during systole is a function of pressure in the opposite ventricle. For the end-systolic relation, the pressure at the end of systole () of both ventricles is formulated as

$$P\_{esLV} = \frac{E\_s \ E\_{lvf}}{E\_s + E\_{lvf}} (V\_{lv} - V\_{l0}) + \frac{E\_{lvf} P\_{rv}}{E\_s + E\_{lvf}},\tag{7.11}$$

$$P\_{es\,RV} = \frac{E\_s \, E\_{rvf}}{E\_s + E\_{rvf}} (V\_{rv} - V\_{r0}) + \frac{E\_{rvf} P\_{lv}}{E\_s + E\_{rvf}}.\tag{7.12}$$

Fig. 7.2: Circuit used to simulate cardiovascular system with interaction from Santamore *et al.*[11]

We then proceeded to extend our model to account for the presence of the pericardium, a fbrous sac that encloses the four chambers of the heart. We followed the implementation from Burkof *et al.* [12] in which the pericardial pressure ( ) is assumed to be equal to an exponential function of the sum of instantaneous LV and RV volumes: = exp[ ( + )].

### **7.2.3 Right Ventricular Remodelling and Barorefex**

The heart continuously regulates its output to match the needs of the body, but it also adapts its anatomy and behavior to changes in pulmonary arterial pressure. During PAH, the heart undergoes long-term deformations in the form of an increase in contractility, and subsequent augment in right ventricular wall volume, . Altered contractility, which is associated with changes in sarcomere length from the ideal optimal value, is evaluated in our model using the formulation from Arts *et al.* [16]. The sarcomere length can be calculated as

$$
\lambda = \left( 1 + 3 \frac{V\_{rv}}{V\_{wall}} \right)^{1/3},
\tag{7.13}
$$

and the change in contractility required to reach such sarcomere length [17] is computed as

$$C = \frac{1+a}{1+a\,e^{b(L\_0 - L\_{s,max})}},\tag{7.14}$$

where = 0.2, = 4−<sup>1</sup> [16], <sup>0</sup> is the desired sarcomere length, and , is the maximum sarcomere length throughout the cycle. Once has been calculated, wall volume () and right ventricular contractility (, ) are respectively updated in an iterative process. Mathematically, a linear approximation from the previous iteration is used in both cases: , = −<sup>1</sup> , and <sup>=</sup> −<sup>1</sup> .

A more short-term cardiac adaptation to changes in heart rate and contractility is barorefex control. This homeostatic mechanism increases the heart rate when arterial pressure decreases, and vice versa [18]. Changes in the heart rate were simulated using the linear approximation [13] = ( −120) +0.855, where is the cardiac period, is the systolic aortic pressure and (the barorefex gain) is fxed at 0.005 /.

### **7.2.4 Atrial Septostomy**

Atrial septostomy is a palliative treatment for medically refractory PAH. The procedure involves using a catheter to intentionally create a hole, called an atrial septostomy defect (ASD), in the interatrial septum. The direct blood fow from the right atrium into the left atrium increases cardiac output, thus alleviating right atrium pressure. As a result, the overall outfow resistance and afterload experienced by the right-side of the heart is reduced, enhancing survival. We account for atrial septostomy in our model by introducing a new lump-parameter for septostomy resistance, , and adding blood fow through the interatrial septal orifce governed by = √ |Δ|, where = <sup>2</sup> is the area of the septostomy in <sup>2</sup> , Δ is the pressure gradient across the orifce in mmHg, and = 2.66.

### **7.2.5 Right Ventricular Assistive Device (RVAD)**

PAH often leads to right ventricular failure, so often assistive devices (RVADs) are required. These devices can be used in a continuous fow capacity, where the fow continuously provides support to the RV independent of time; or as a pulsating RVAD, where the fow is time dependent and in sync with the cardiac cycle. Both continuous and pulsatile fow pumps have previously been tested computationally and experimentally [9, 19]. Additionally, the RVAD can be implemented in either the right ventricle or right atrium. To implement the RVAD into the model, we add an additional fow, . The fow into the pulmonary artery will then be given by

$$\frac{dV\_{pa}}{dt} = \mathcal{Q}\_{pa} - \mathcal{Q}\_{pv} + \mathcal{Q}\_{rvad}.\tag{7.15}$$

RVAD will take fuid from the right ventricle or atrium, which results in either

$$\frac{dV\_{rv}}{dt} = \mathcal{Q}\_{tr} - \mathcal{Q}\_{pvv} - \mathcal{Q}\_{rvad} \tag{7.16}$$

or

$$\frac{dV\_{ra}}{dt} = \mathcal{Q}\_{sv} - \mathcal{Q}\_{tr} - \mathcal{Q}\_{rad},\tag{7.17}$$

respectively. The form of depends on the type of RVAD. The continuous fow RVAD only depends on the diference in pressure between the source (the right atrium or ventricle) and the pulmonary artery. Analytically, this can be approximated using a linear relationship [8], where is time independent. We also implement the pulsatile fow pump as described in Gohean *et al.* [9], which has demonstrated better performance compared to the continuous fow pumps [19], when properly synchronized. Now is time dependent and follows the formulation in Gohean *et al.* [9], given as

$$Q\_{r\text{vad}}(t,\Delta P) = \begin{cases} \frac{2V\_s}{T\_l} \left(\frac{1}{2} - \frac{1}{2}\cos\left(\frac{2\pi(t - T\_d)}{T\_l}\right)\right) & \text{if } t \in [T\_d, T\_d + T\_l],\\ 0 & \text{otherwise.} \end{cases} \tag{7.18}$$

where is the time delay from the start of the cardiac cycle, is the duration of the pulse, and is the stroke volume in . The time, , is normalized to the start of the cardiac cycle.

### **7.3 Results**

We successfully implemented a ventricular hemodynamic model and simulated PAH based on data from Punnoose *et al.* [8] to study the efects of this aggressive disease on RV hemodynamics. The right and left ventricular hemodynamic parameters reproduced by our model are reported in Table 7.1.

Table 7.1: Stroke volume (SV), cardiac output (CO), and ejection fraction (EF) from the model for healthy, severe PAH, continuous RVAD treatment, and pulsatile RVAD treatment. Stroke volume is given in and cardiac output in /.


An important indicator of the severity of PAH is RV function and morphology [20]. Figure 7.3 shows pressure-volume loops (PV loops) corresponding to various severities of PAH. The shapes of the PV loops, which provide information on cardiac performance, load and coupling, change with progression of PAH. To simulate the diferent stages of PAH, the elastance, resistance and compliance parameters as well as the heart rate were modifed based on values from previous studies [8]. As the disease progresses, resistance in the pulmonary artery increases, heart rate increases, and end-systolic elastance increases in the left ventricle, while it decreases in the right ventricle. For the LV, the normal shape of the PV loop is rectangular, whereas it is usually rounded for the RV. For both ventricles, the PV loops become more narrow, which is indicative of reduced stroke volume.

Fig. 7.3: PV loops for healthy and diseased hearts, from mild PAH to PAH-induced heart failure. Dark blue represents the most severe case, cardiogenic shock (CGS).

### **7.3.1 Sensitivity Analysis**

A sensitivity analysis was a useful tool to determine the correlation between input model parameters and the resulting hemodynamic features of the system. By altering each model parameter, one at a time, we were able to better interpret the model outputs when the inputs were changed in PAH treatment simulations. The model behavior changed drastically when the pulmonary arterial resistance () was altered, as illustrated in Figure 7.4. The values plotted correspond to 0.25, 0.50, 0.75 and 1 mmHg.s/mL. As increased (as in PAH), the stroke volume decreased substantially. The total mechanical energy generated by ventricular contraction also decreased, since the area encompassed within the PV loops (i.e., the stroke work) decreased with higher values. In addition to pulmonary arterial resistance, other parameters that also had a considerably high impact on the model outputs were RV end-systolic elastance and the exponential parameter () for the right ventricular end-diastolic PVR ( 7.5).

Fig. 7.4: PV loops for several degrees of pulmonary arterial resistance.

### **7.3.2 Ventricular Interaction and the Pericardium**

As explained in subsection 7.2.2, we aimed to improve the realism of the model by including the interaction via the interventricular septum. We implemented three timevarying elastances with the same values used by Santamore *et al.* [11], and simulated diferent degrees of PAH by changing the values of resistance as in Punnoose *et al.* [8], and keeping the same values of elastance from Santamore *et al.* [11]. Figure 7.5 shows the diferent pressure-volume (PV) loops of the left and right ventricles when the interventricular interaction is turned on. The overall behavior of the PV loops is the same as those without interaction: there is a leftward shift of the LV PV loop and a rightward shift of the RV PV loop with increasing PAH severity. This is likely due to the RV expanding and pushing against the LV, causing the right ventricular volume to increase and the left ventricular space to decrease, eventually leading to ventricular septum bowing. Compared to Figure 7.3, Figure 7.5 displays PV loops that shift more uniformly along an axis as severity increases. One possible explanation for this behaviour is the use of the same elastance values for the ventricular interaction simulations at each severity level. This behaviour, however, is not as strongly evident as in the previous simulations, which could be explained by the fact that we kept the same elastance from Santamore *et al.* [11] for each PAH degree, and because we used a combination of two parameter sets (Punnoose *et al.*[8] and Santamore [11]). We can still conclude that the introduction of the interventricular interaction does not lead to strong changes for this model. This result is in agreement with other studies that report small efects from the septum [11].

Fig. 7.5: Pressure-Volume loops with diferent degrees of PAH when the ventricular interaction via the septum is included. Dark blue represents cardiogenic shock (CGS).

Secondly, the presence of the pericardium represents an additional constraint to both ventricles. Figure 7.6 shows the efect of the pericardium in a healthy PV loop. When this element was added to our model, it lead to a reduction in the overall left ventricular stroke volume and an increase in the end-diastolic pressure of the right ventricle. The PV loop shifts slightly upward of the end-diastolic pressure-volume relationship (EDPVR), which is representative of the heart's ventricular passive compliance curve. In addition, the stroke volume is smaller and the peak pressure is lower in presence of pericardium. The stroke volume reduction observed in our model is consistent with physiological and clinical observations that the pericardium constricts heart dilation in the diastolic phase.

Fig. 7.6: Pressure-Volume loops of healthy condition with and without pericardium

### **7.3.3 Remodelling the Right Ventricle**

Including the remodelling efects of the RV resulted in improved performance of the heart. Computing the change in contractility, as described by Equation 7.14, yielded improved model output for an increase in pulmonary pressure (Figure 7.7). Additionally, the arterial baroreceptor refex system changed the model behavior in presence of PAH. The heart rate and contractility (heart rate efect shown in Figure 7.8), which are two parameters that have a signifcant efect on the model's output, are adjusted by the barorefex to control fuctuations in aortic blood pressure. Thus, impaired barorefex control of the heart in PAH would typically contribute to worse outcome of the disease.

Fig. 7.7: Efect of remodelling the RV on the pressure volume relationship in the left and right ventricle.

Fig. 7.8: Heart rate control, characteristic of barorefex, has a high impact on the shapes of the PV loops.

### **7.3.4 Atrial Septostomy**

Consistent with the left- and right-ventricular hemodynamics behavior observed in Punnoose *et al.*, our atrial septostomy defect simulation with resulting right-to-left shunt yields minimal leftward shifts in both the right atrial and right ventricular pressure-volume (PV) loops toward lower volumes and pressures [8]. Simultaneously, a rightward shift in the left atrial and left ventricular PV loops refect increased flling pressures and end-diastolic volumes as shown in Figure 7.9.

Based on the results from our model, atrial septostomy signifcantly increased left ventricular stroke volume as a result of elevated end-systolic pressure and expanded end-diastolic volume. Meanwhile, the left ventricular end-diastolic pressure and end-systolic volume experienced no signifcant change with atrial septostomy. Collectively, these results indicate atrial septostomy causes an increase in left ventricular cardiac output with negligible change in the systemic vascular resistance. The right ventricular pressure, volume, and stroke volume remained largely the same with or without atrial septostomy and are thus determined to be minimally responsive to atrial septostomy treatment.

Prior to atrial septostomy treatment, the left and right atria are isolated, pressurized elastic chambers, with the right atrium having a higher pressure than the left atrium. The chamber pressure values are within reasonable ranges of 1 to 7 mmHg for the right atrium and 1 to 6 mmHg for the left atrium, as the atria are low pressure chambers meant for gathering blood from veins. After atrial septostomy treatment, however, we observed an increase of the end-systolic pressure and the end-diastolic volume in the left atrium as well as a slight decrease in the right atrium end-systolic pressure and end-diastolic volume. This is likely due to the two atria being connected via the atrial septal defect post-surgery, driving the pressures in the left and right atria to an equilibrium.

Our model observed a decrease in interatrial pressure gradient. This is demonstrated in Table 7.2 by a 34% elevation in left atrial end-diastolic pressure to match an 11% decrease in right atrial end-diastolic pressure. This outcome aligns with our prediction that the two atrial chambers will reach a pressure equilibrium when two heart chambers are connected via the atrial septal defect.

Fig. 7.9: Efect of atrial septostomy (ASD) on the pressure-volume relationship in each heart chamber.



### **7.3.5 Right Ventricular Assistive Device**

We fnd that our RVAD implementation behaves similarly to the one in Punnoose *et al.* [8] and Gohean *et al.* [9] (Table 7.1 and Figure 7.10). Simulated PAH reduced left ventricular cardiac output, stroke volume and ejection fraction (Table 7.1). This is evident in Figure 7.10, with a leftward shift in the left ventricular PV loops indicating left ventricular unloading. This is improved by pulsatile and continuous RVAD at a fow rate of 3.5 L/min from the right atrium or ventricle, indicated by a rightward shift in the left ventricular PV loops for both treatments, with the pulsatile fow performing better than the continuous fow RVAD (Figure 7.10 right). Conversely, the right ventricular end-systolic and diastolic pressures increase in the simulated PAH. As shown in Figure 7.10 (left), a rightward shift in the PAH PV loop is seen, which is not improved by either the pulsatile or continuous RVAD.

Fig. 7.10: Ventricular PV loops in healthy, severe PAH and treatment cases. The right ventricular PV loops for severe PAH shift to the left, but are partially recovered by pulsatile or continuous RVAD (left). Left ventricular PV loops in healthy, PAH and treated simulations of hemodynamics (right).

### **7.4 Discussion**

In this study, we have modeled the infuence of PAH on the mechanics of the circulatory system at various stages of the disease, based on the results from Punnoose *et al.*[8]. A sensitivity analysis demonstrated the signifcant dependence of our model on the resistance in the pulmonary artery, but many factors beyond resistance, such as contractility and stifness, have an efect on the model too. To improve upon the base model, the addition of the interaction between the ventricular walls through the septum [11] and the efects of the pericardium [13] were included. We fnd a larger infuence from the pericardium, which agrees with previous results [11, 13]. Further, we investigated long and short term adaptions of the heart to the increase in pressure of the pulmonary artery, as well as the efect of treatments to PAH. Specifcally, the target of our computational study was to provide insight into how atrial septostomy and RVAD afect heart function [8, 9].

Our model exhibited hemodynamic efects of atrial septostomy similar to those observed in Punnoose *et al.* for all four heart chambers. This validated our implementation to study the efect of atrial septostomy on PAH-induced right-side heart failure [8]. Results from our simulation indicate that atrial septostomy indeed restores the heart to a healthier functional state under the assumption of *ceteris paribus*, consistent with clinical observations. The end-diastolic pressure-volume relationship (EDPVR) on the PV loop represents the passive compliance of the heart ventricles. With atrial septostomy, the EDPVR remained the same in both the left and right ventricles, implying that the heart ventricular tissue stifness is not afected by the treatment. The end-systolic pressure-volume relationship (ESPVR), which represents ventricular tissue contractility, remained the same for the left ventricle and decreased for the right ventricle with atrial septostomy treatment, denoting less contractility. The introduction of atrial septostomy in the interatrial septum likely decreased the right ventricle's need to contract to maintain the same cardiac output under PAH. Collectively, these data suggest that, while atrial septostomy does not change the myocardium intrinsic biomechanical properties, the surgical procedure seems to promote cardiac output, improve heart performance, and reduce heart failure symptoms. Summed together, our model provided results consistent with current clinical understanding of atrial septostomy efects on cardiovascular biomechanics, physiology, and literature.

The RVAD resulted in improved cardiac performance for a severe PAH setting by restoring cardiac output and ejection fraction in the LV. In our simulations, we observe better performance from the pulsatile fow compared to the continuous fow pumps. Despite the improved performance, the pulsatile fow pumps typically sufer from durability issues [9], so continuous pumps must be investigated too. Using this model, the infuence of both types of pumps on the hemodynamics can continue to be studied in future research.

Our model was able to capture the efects of PAH on heart function, but there are some limitations to our work. Our ability to investigate the growth and remodelling was limited by the fact that we didn't account for the three-dimensional structural efects from the heart chambers, which alter the anatomy and behavior of the RV. To overcome this limitation, our framework could incorporate a 3D fnite element model of the ventricles [21]. This extension of the lumped-parameter model would provide a more detailed characterization of ventricular mechanics as well as a more realistic representation of ventricular interaction via the septum. On the other hand, the limited literature sources with specifc experimental data relevant to our study partly afect the accuracy of our results. There is signifcant variability in the biomechanical data between diferent cardiac studies, which makes it challenging to incorporate elements from multiple studies and unify them in one model [8, 11]. Conducting our own *in vivo* experiments could be useful to obtain more reliable parameters for the model. That is, experimental measurements would be a key next step to corroborate the results and develop an improved version of our model.

### **7.5 Conclusion**

In this computational study, we have implemented a closed-loop lumped-parameter model of the heart and circulatory system for studying the cardiovascular adaptations in the RV to PAH. We demonstrated that RVAD and atrial septostomy are efective treatments to improve left ventricular flling and cardiac output. We also extended our base model to account for direct ventricular interaction via the interventricular septum and the presence of the pericardium, showing that the pericardium has a larger efect relative to the ventricular interaction. Our present implementation lays the foundation for the addition of more components to the model that would augment its realism even further. We believe that highly complex computational models will continue to be critical to achieve our vision of dynamically simulating and characterizing ventricular mechanics in PAH.

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.